AEROSPACE  REPORT  NO. 
TR-2010(3905)-2 


Trapped  Electron  Model  2  (TEM-2) 


April  25.2010 


Paul  O'Brien  and  Timothy  B.  Guild 
Space  Science  Applications  Laboratory 
Physical  Sciences  Laboratories 


Prepared  for: 

National  Reconnaissance  Office 
14675  Lee  Rd. 

Chantilly,  VA  20151-1715 


Authorized  by:  National  Systems  Group 


APPROVED  FOR  PUBLIC  RELEASE: 
DISTRIBUTION  UNLIMITED 


This  report  was  submitted  by  The  Aerospace  Corporation,  El  Segundo,  CA  90245-4691 , 
under  Contract  No.  FA8802-09-C-0001  with  the  Space  and  Missile  Systems  Center,  483  N. 
Aviation  Blvd.,  El  Segundo,  CA  90245.  It  was  reviewed  and  approved  for  The  Aerospace 
Corporation  by  J.  H.  Clemmons,  Principal  Director,  Space  Science  Applications  Laboratory; 
and  G.  A.  Davis,  General  Manager,  Imagery  Programs  Division.  David  Byers  was  the 
project  officer  for  the  program. 

This  technical  report  has  been  reviewed  and  is  approved  for  publication.  Publication  of  this 
report  does  not  constitute  NRO  approval  of  the  report's  findings  or  conclusions.  It  is  pub¬ 
lished  only  for  the  exchange  and  stimulation  of  ideas. 


David  Byers 


IMINT  R&T 


SC-1 865(2,  5666,  61,  JS) 


REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  date  sources,  gathenng  end 

maintaining  the  data  needed  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information 
including  suggestions  for  reducing  this  burden  to  Department  of  Defense,  Washington  Headquarters  Services  Directorate  for  Information  Operations  end  Reports  (0704-0188)  1215  Jefferson 

Davis  Highway  Suite  1204,  Arlington,  VA  22202-4302  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  lew  no jierson  shall  be  subject  to  eny  penalty  for  failing  to 
comply  with  a  collection  of  information  if  it  does  not  display  a  currently  valid  OMB  control  number  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 

1 .  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

25-04-2010 

3.  DATES  COVERED  (From  -  To) 

4.  TITLE  AND  SUBTITLE 

Trapped  Electron  Model  2  (TEM-2) 

5a.  CONTRACT  NUMBER 

FA8802-09-C-000 1 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

Paul  O’Brien  and  Timothy  B.  Guild 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

T  he  Aerospace  Corporation 

Physical  Sciences  Laboratories 

El  Segundo,  CA  90245-4691 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

TR-2010(3905)-2 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

National  Reconnaissance  Office 

14675  Lee  Rd. 

Chantilly,  VA  20151-1715 

10.  SPONSOR/MONITOR’S  ACRONYM(S) 

NRO 

11.  SPONSOR/MONITOR’S  REPORT 

NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited. 

13.  SUPPLEMENTARY  NOTES 


2.0|  20  <420  MS' 

14.  ABSTRACT 

We  present  a  next-generation  radiation  specification  for  the  electron  radiation  belts.  The  specification  includes  an  empirical 
model  of  the  statistical  variation  and  spatiotemporal  cov  ariation  of  electron  fluxes  from  40  keV  to  10  MeV  in  the  inner  and  outer 
zones,  as  well  as  a  dynamic  full  solar-cycle  reanalysis,  which  provides  a  global  snapshot  of  the  electron  fluxes  every  6  hours 
from  1992  to  2008.  The  statistical  model  is  derived  from  data  taken  by  the  S3-3,  CRRES,  Polar,  and  SCATHA  spacecraft.  The 
reanalysis  is  driven  with  data  from  HE:0-1 ,  HEO-3,  ICO,  and  SAMPEX.  We  validate  the  reanalysis  by  comparing  a  test  reanal¬ 
ysis  without  ICO  data  to  actual  ICO  observations.  T  he  statistical  model  alone  can  demonstrate  Monte  Carlo  characterization  of 
uncertainties  in  radiation  specifications  for  internal  charging  and  total  ionizing  dose  due  to  electrons.  The  reanalysis  can  demon¬ 
strate  engineering  applications,  such  as  a  “standard  solar  cycle”  for  spacecraft  design,  as  well  as  scientific  applications,  such  as 
solar  wind  correlation  studies,  initial  and  boundary  conditions  for  numerical  simulations,  and  principal  component  analysis.  We 
demonstrate  application  of  the  model  to  radiation  specifications  and  to  principal  component  analysis  of  the  electron  belts.  Due  to 
the  prototype  nature  of  the  model,  it  is  insufficiently  accurate  for  definitive  conclusions  to  quantitative  investigations. 

15.  SUBJECT  TERMS 

Space  environment,  Radiation  belts.  Radiation  specification  models,  Radiation  belt  climatology.  Reanalysis 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

18.  NUMBER 

19a.  NAME  OF  RESPONSIBLE 

OF  ABSTRACT 

OF  PAGES 

PERSON 

Paul  O’Brien 

a.  REPORT 

UNCLASSIFIED 

b.  ABSTRACT 

UNCLASSIFIED 

C.  THIS  PAGE 

UNCLASSIFIED 

Leave  blank 

61 

19b.  TELEPHONE  NUMBER 
(include  area  code) 

571-307-3978 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  239.18 


Ackmowledgments 


This  work  was  supported  under  The  Aerospace  Corporation's  Independent  Research 
and  Development  Program,  by  NASA  TWS  Targeted  Research  and  Technology 
Grant  NNG05GM22G,  and  by  the  NRO’s  Proton  Spectrometer  Melt  Research 
program.  A  merged  Polar/CEPPAD  dataset  was  provided  by  R.  Friedel  (LANL).  An 
L-sort  version  of  the  SAMPEiX  PET/ELO  data  was  provided  by  S.  Kanekal  (LASP). 
A  revised,  temperature-corrected  version  of  the  CRRES/EIEEF  data  was  provided  by 
G.  Ginet  (AFRI>  and  MIT/Eincoln  Tab).  We  made  extensive  use  of  the  1RBELM 
(formerly  ONERA)  library  of  magnetic  field  models  and  field  line  tracing  routines. 
The  ICO  data  was  funded  by  AFRL,  and  the  EIEO  and  SCATHA  data  were  provided 
by  our  colleagues  at  The  Aerospace  Corporation.  The  authors  thank  M.  Eooper  for 
providing  calculated  response  curves  for  the  HEO  and  ICO  channels.  We  thank  A. 
Vampola  for  follow-up  discussion  regarding  his  S3-3,  CRREiS,  and  other  magnetic 
spectrometers.  In  addition,  the  authors  would  like  to  thank  J.  Fennell,  J.  Roeder,  J.  B. 
Blake,  J.  Mazur,  S.  Huston,  S.  Bourdarie,  D.  Boscher,  R.  Selesnick,  D.  E*rautigam,  C. 
Roth,  R.  Quinn,  P.  Whelan,  W.  R.  Johnston,  C.  Lind  Strom,  J.  Albert,  K.  Perry,  B. 
Wie,  D.  Byers,  R.  Weigel,  I).  Kondrashov,  Y.  Shprits,  D.  Vassiliadis,  J.  Roller,  and 
G.  Reeves  for  helpful  discussions. 


ill 


Contents 


1  Introduction  1 

1.1  Energetic  electron  hazards  to  spacecraft .  1 

1.2  Review  of  established  specification  models .  1 

1.3  Review  of  other  next  generation  specification  models .  2 

1.4  A  different  approach:  A  statistical  model  .  3 

2  Statistical  prior  model  4 

2.1  Data  used  to  construct  the  statistical  prior  model .  8 

2.2  Data  processing  to  develop  the  prior  model .  9 

2.2.1  Binning .  9 

2.2.2  Daily  averaging  within  bins .  11 

2.2.3  Computing  9  and  c,o\  80  in  each  bin .  13 

2.2.4  Interpolating  to  the  model  grid .  13 

2.2.5  Flux  map  diagnostics  .  15 

2.2.6  Computing  the  errors  in  the  interpolated  flux  map  0 .  18 

2.2.7  Computing  spatial  and  spatiotemporal  covariances .  19 

2.2.8  Principal  components . 21 

2.3  Engineering  application:  Monte  Carlo  mission  scenarios . 22 

2.3.1  GEO  example . 23 

2.3.2  GPS  example . 25 

2.3.3  Dose  depth  curves . 28 

3  Reanalysis  procedures  30 

3.1  Stating  the  data  assimilation  problem . 30 

3.2  The  likelihood  to  be  maximized . 31 

3.3  The  statistical  prior  model  as  a  likelihood  function . 31 


IV 


3.4  1  he  grid  interpolation  functions . 33 

3.5  Measurement  functions . 34 

3.6  Pseudo- Poisson  error  distribution . 34 

3.7  Gaussian  error  distribution  . 35 

4  Reanalysis  data  36 

4.1  11  EC)  and  ICO  data . 36 

4.2  SAMPEX  data . 38 

4.3  Accounting  for  latent  data . 38 

5  Reanalysis  results  43 

5.1  Validation  and  known  issues . 43 

5.2  Engineering  application . 46 

5.3  Science  application . 46 

5.3.1  Empirical  determination  of  the  time  evolution  operator . 46 

5.3.2  Posterior  principal  component  analysis . 47 

6  Summary  56 

A  TEM-1  57 

B  Useful  derivatives  58 


v 


Figures 


1  An  example  of  the  Weibull  distribution  fit  to  observed  fluxes  in  a  selected  energy  channel 
in  a  bin  of  aeq  and  Lm.  The  distribution  fits  well  everywhere  except  the  very  lowest  2-3% 
of  the  distribution,  which  is  of  little  consequence  for  spacecraft  effects.  The  top  panel  is  a 
quantile-quantile  plot,  and  the  bottom  panel  shows  the  tails  of  the  cumulative  distribution 

and  exceedance  probability .  6 

2  A  pairwise  scatter  plot  of  the  transformed  fluxes,  u,  for  two  bins,  superimposed  on  the 
probability  density  of  the  Gaussian  copnila  with  correlation  parameter  0.71  taken  from  the 
data.  Contour  lines  are  provided  at  copula  densities  of  0  to  4  in  steps  of  0.5.  As  expected, 


the  points  cluster  where  the  density  is  high  (red) .  7 

3  Flow  chart  depicting  data  processing  steps  that  lead  to  the  prior  model  data  tables.  ...  10 

4  TEM-2  data  coverage  combined  over  all  missions  in  selected  energy  ranges . 12 


5  Median  and  95th  percentile  flux  maps  for  one  OR  RES  MEA  channel,  along  with  errors 
on  6  and  coverage.  A  black  “X”  indicates  that  an  aeq  —  Lm  bin  did  not  meet  the  bin 


selection  criteria  defined  in  section  2.2.1,  usually  because  of  too  few  days  of  data .  14 

6  TEM-2  maps  of  median  and  95th  percentile  fluxes .  16 

7  Comparison  of  TEM-2  to  AES  and  CRRESELE  along  an  equatorial  radial  profile . 17 

8  Comparison  of  TEM-2  to  AES  and  CRRESELE  along  a  field  line  through  (4. 5,0,0)  R&  in 

the  outer  zone .  17 

9  Comparison  of  TEM-2  to  AES  and  CRRESELE  along  a  LEO  t  rajectory  through  the  South 

Atlantic  Anomaly .  IS 


10  Median  (top)  and  95th  percentile  (bottom)  fluxes.  Data  from  in  situ  observations  in  the 
bin  are  plotted  in  color.  The  nearest  neighbors  fit  is  indicated  in  black.  Also  shown  is  the 
“base  0 ”  given  by  22.  Error  bars  on  the  points  and  the  fit  indicate  one  standard  deviation.  20 

11  Comparison  of  TEM-2  scenarios  and  reanalysis  to  AES  and  GOES-10  observations  for  a 

real  interval.  GOES- 10  flux  is  multiplied  by  37r  to  approximate  conversion  from  unidi¬ 
rectional  to  omnidirectional  flux,  assuming  that  the  flux  follows  a  sin'J/ 4  a  distribution 
in  local  pitch  angle.  CRRESELE,  as  plotted,  is  driven  by  uApl5”,  the  15-dav  trailing 
average  of  the  Ap  index  [Brautigam  and  Bell ,  1995].  The  percentiles  shown  are  taken 
from  200  Monte  Carlo  scenarios,  of  which  only  20  are  shown  (gray) . 23 

12  Comparison  of  TEM-2  scenarios  and  GOES-10  observed  statistical  distribution  of  >  2 
MeV  electron  flux.  The  plot  format  follows  that  of  Figure  1.  GOES-10  flux  is  multiplied 
by  37r  as  in  Figure  11.  The  TEM-2  distribution  shown  is  taken  from  200  Monte  Carlo 


scenarios . 24 

13  10-year  environments  at  GEO.  Top:  Mean  flux  from  various  models;  AES  (MIN  and  MAX 

are  the  same  at  GEO).  Bottom:  worst-case  24-hour  exponentially- weighted  flux  average.  .  26 


vi 


14  10-year  environments  at  GPS.  Top:  Mean  flux  from  various  models;  AE8.  Bottom:  worst- 

case  24-hour  exponentially- weighted  flux  average . 27 

15  10-year  total  dose  environments  at  GEO  (top)  and  GPS  (bottom) . 29 

16  Calculated  energy  response  curves  for  HEO  and  ICO  channels . 37 

17  (Top)  SAMPEX  PET/ELO  intensity  can  be  used  as  a  proxy  for  HEO-3  E5.  (Bottom) 

Pile  parameters  of  the  proxy  relationship  vary  with  L  shell . 39 

18  Autocorrelation  for  HEO-3  electron  fluxes  at  1  day  lag  for  all  orbit  status  and  energv 

channels . 41 

19  Autocorrelation  for  HEO-3  electron  fluxes  at  2  days  lag  versus  the  square  of  the  correlation 
at  1  day  lag  for  all  orbit  status  and  energy  channels.  The  good  agreement  suggests  that 


the  HEO-3  autocorrelation  is  a  first-order  autoregressive  process . 42 

20  Flux  versus  Lm  profiles  versus  time  in  the  FEM-2  Reanalysis:  (top)  the  entire  reanalysis, 

(middle)  one  year,  centered  on  the  July  2004  storm,  and  (bottom)  one  month  centered  on 
July  24,  2004. . .  44 

21  Validation  of  TEM-2  reanalysis  versus  ICO  (in  and  out  of  sample) . 45 

22  The  first  posterior  principal  component  of  t  lie  TEM-2  reanalysis,  for  equatorially  mirroring 

electrons . 48 

23  the  second  posterior  principal  component  of  the  1  EM-2  reanalysis.  Top:  equatorially 

mirroring  electrons.  Bottom:  selected  isocontours . 50 

21  I  he  third  and  fourt  h  posterior  principal  components  of  the  FEM-2  reanalysis,  for  equa¬ 
torially  mirroring  electrons . 51 

25  Top:  TEM-2  reanalysis  posterior  principal  components  during  the  October  1998  GEM 


storm.  Lower  panels:  observed  electron  flux  at  GEO,  observed  solar  wind  speed,  and  the 


magnetic  storm-time  disturbance  index  I)st . 52 

26  1  EM-2  reanalysis  posterior  principal  components  and  measurements  during  the  October 

2000  GEM  storm,  in  the  format  of  Figure  25 . 53 

27  Prediction  efficiency  of  1  EM-2  reanalysis  posterior  principal  components  with  respect  to 

GOES  >  2  MeV  electron  flux,  Dst1  and  solar  wind  speed  (Vste) . 51 


28 


Impulse  response  of  TEM-2  reanalysis  posterior  principal  components  with  solar  wind 
speed . 


Tables 


1  Data  used  to  build  the  statistical  prior  model .  9 

2  Coefficients  of  0^°^  .  14 

3  Data  channels  used  to  drive  the  reanalysis . 36 

4  Coefficients  for  1-day  autocorrelation  function  vs  L  for  I1EO-3  fluxes . 40 

5  ICO  out-of-sample  validation  statistics . 43 


viii 


1  Introduction 


We  have  developed  a  prototype  next  generation  radiation  specification  for  electrons  with  energies  from 
10  keV  to  10  MeV  trapped  in  the  Earth’s  magnetic  field.  We  denote  this  model  Trapped  Electron 
Model  2  (TEM-2).  Like  its  immediate  predecessor,  TEM-1  (unpublished,  see  Appendix  A),  TEM-2  was 
built  as  a  prototype  for  future  design  specifications  models.  However,  in  addition  to  its  value  as  a 
prototype,  TEM-2  also  facilitates  some  new  kinds  of  scientific  research  in  the  electron  belts.  We  will 
demonstrate  both  engineering  and  scientific  applications  of  TEM-2.  TEM-2  is  an  alpha  version  of  the 
upcoming  AE9  model.  Although  it  demonstrates  new  scientific  and  technical  capabilities,  it 
incorporates  some  known  but  unresolved  issues  with  the  underlying  data  (see  section  2.1);  therefore 
quantitative  outputs  of  TEM-2  are  too  preliminary  to  be  used  for  satellite  design  or  quantitative 
science. 


We  will  begin  with  a  review  of  the  importance  of  energetic  electrons  for  spacecraft  design  and 
operations,  and  then  we  will  review  the  existing  suite  of  specification  models,  including  other  next 
generation  models  which  are  in  development.  Next,  we  will  describe  our  own  model  and  its 
development  in  detail.  We  will  also  demonstrate  the  model  in  an  engineering  application  and  in  a 
scientific  one. 


1.1  Energetic  electron  hazards  to  spacecraft 

Energetic  electrons  trapped  in  the  Earth’s  radiation  belt  typically  pose  two  hazards  to  spacecraft:  total 
ionizing  dose  and  internal  charging.  The  total  ionizing  dose  hazard  arises  from  the  degradation  of 
electronics  and  materials  from  accumulated  absorption  of  ionizing  radiation.  Spacecraft  designers  must 
consider  ionizing  dose  from  keV  to  MeV  electrons  when  choosing  parts  and  shielding  to  reduce  the 
likelihood  of  premature  failure  of  a  spacecraft  or  subsystem  [e.g.,  Tribble  et  a/.,  1999].  Internal  charging 
by  energetic  electrons  is  one  potential  cause  of  electrostatic  discharge  (ESI)).  Energetic  electrons  can 
penetrate  spacecraft  shielding  and  deposit  their  charge  internal  to  the  spacecraft  surface.  If  a  large 
electrostatic  potential  builds  up  before  the  charge  can  dissipate,  an  ESI)  may  occur,  leading  to 
phantom  commands,  short  circuits,  and  a  host  of  other  kinds  of  anomalies  [e.g.,  NASA-HDBK-^002, 
1999].  "Therefore,  a  spacecraft  designer  needs  a  specification  of  the  mean  electron  environment  which 
leads  to  total  dose  specifications,  and  the  worst  case  environment,  which  leads  to  internal  charging 
specifications.  Because  the  mission  mean  and  worst  case  environments  carry  uncertainty  due  to 
radiation  belt  dynamics  and  measurement  errors,  there  is  need  for  next  generation  radiation 
specification  models  that  can  describe  the  dynamic  variability  of  the  environment  and  characterize 
errors  in  the  models  themselves  [e.g.,  Vette,  1978;  Gussenhoven  et  al. ,  1994;  Daly  et  ai,  1996]. 

1.2  Review  of  established  specification  models 

The  industry  standard  for  electron  radiation  specification  is  the  NASA  AES  model  [Vette,  1991].  This 
model  provides  global  coverage  of  the  inner  and  outer  belt  electrons  (Mcllwain  Lm  from  1.2-11)  with 
energies  from  0.04  to  7  MeV.  Rudimentary  dynamic  variation  in  the  environment  is  provided  in  the 
form  of  different  specifications  for  solar  maximum  and  solar  minimum. 

The  AFRL  CRRESELE  model  [ Brautigam  and  Bell ,  1995]  describes  electrons  with  energies  from  0.5  to 
fi.fi  MeV,  Lm  2. 5-6.8,  and  magnetic  latitude  below  68  degrees.  CRRESELE  provides  some  indication  of 
the  dynamic  environment  in  the  form  of  both  average  and  worst-case  radiation  maps,  as  well  as 
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distinct  radiation  maps  for  various  bins  of  the  15-day  average  of  the  Ap  magnetic  index.  Presumably 
because  of  its  limited  spatial  and  energy  coverage,  CRRESELE  has  not  been  widely  adopted. 

Finally,  Vampola  [1996]  produced  an  update  to  AE8-MIN  outer  zone  flux  map  using  data  from  the 
MEA  instrument  on  CRRES.  Vampola  developed  a  neural  network  that  specified  the  outer  zone 
electron  fluxes  based  on  the  preceding  few  daily  averages  of  the  Kp  index.  He  then  ran  this  neural 
network  model  over  the  entire  multi-decade  history  of  Kp.  He  selected  solar  minimum  years  from  the 
results  and  computed  an  average  outer  zone. 

There  is  an  extensive  literature  identifying  various  shortcomings  of  these  models  [e.g.,  Vette,  1978; 
Gussenhoven  et  ai,  1994;  Daly  et  a/.,  1996].  In  this  paper,  we  are  concerned  mainly  with  the  lack  of  a 
model  with  sufficient  energy  and  spatial  coverage  that  can  also  provide  estimates  of  the  dynamic 
variability  of  the  environment  and  errors  in  the  model  itself.  We  note  that  at  the  time  AE8  and  the 
CRRES  models  were  built,  there  were  few  to  no  data  sets  of  sufficient  temporal  and  energy /spatial 
coverage  to  properly  describe  the  global  electron  dynamics. 


1.3  Review  of  other  next  generation  specification  models 


The  development  of  a  Monte  Carlo  scenario  capability  for  a  radiation  belt  model  is  presently  unique  to 
our  efforts.  However,  the  reanalysis  aspect  of  our  effort  is  similar  to  several  others. 

We  are  aware  of  other  efforts  to  develop  reanalysis  models  of  the  electron  radiation  belts.  One  effort,  at 
ONER  A  (Office  National  d?  Etudes  et  Recherches  Aerospatiales)  builds  on  the  well-known  Salammbo 
model  [Maget  et  a/.,  2007].  By  adding  data  assimilation  of  measurements  from  LANL  sensors  on  GPS 
and  geosynchronous  (GEO)  satellites,  the  ONERA  group)  has  obtained  promising  results,  including  a 
first  physics-based  data  assimilation  of  solar-cycle  duration  (a  “reanalysis”).  At  present,  the  data 
assimilation  scheme  being  used  is  direct  insertion,  which  seems  to  work  well  on  multi-day  timescales 
because  diffusive  transport  quickly  smoothes  out  the  sharp  features  produced  by  direct  insertion.  This 
effort  has  produced  a  model  outer  zone  average  on  a  year-by-year  basis  for  a  full  solar  cycle  [ Bourdarie 
et  a/.,  2009]. 

A  second  group,  the  LANL  “DREAM”  team,  is  pursuing  a  more  ambitious  model  that  couples  the 
radiation  belt  into  a  global  model  that  includes  the  ring  current,  plasmasphere  and  convection  electric 
fields  [Koller  et  ai ,  2007].  At  time  of  writing,  DREAM  has  produced  promising  initial  results  from  an 
uncoupled  outer  zone  model  that  employs  an  ensemble  Kalman  filter  for  data  assimilation  of  a  radial 
diffusion-only  simulation.  This  model  is  driven  bv  LANL  GEO  and  GPS  data  and  by  data  from 
NASA’s  Polar  spacecraft. 

A  third  group,  the  UCLA  “VERB”  team,  is  also  developing  a  physics-based  data  assimilation  model. 
The  VERB  model  has  recently  been  used  to  compare  reanalyses  from  CRESS  and  Akebono 
observations  [ Ni  et  ai ,  2009;  Shprits  et  ai ,  2007]. 

We  maintain  frequent  contact  with  the  ONERA,  LANL,  and  UCLA  teams,  so  that  our  efforts  and 
theirs  can  complement  each  other.  We  collaborate  extensively  on  data  intercalibration,  since  that  work 
would  otherwise  be  duplicated  wastefully.  We  feel  that  ultimately  data  assimilation  into  a  physical 
model  (like  theirs)  will  supersede  our  statistical  reanalysis,  so  that  our  results  can  be  seen  as  an  early 
step  in  the  process  and  as  a  potential  resource  for  future  improvement  of  their  physics-based  reanalysis. 
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1.4  A  different  approach:  A  statistical  model 


Unlike  the  models  described  in  section  (1.3),  our  model  does  not  include  a  physics-based  numerical 
simulation.  Typically,  observations  do  not  fully  specify  the  global  state  of  the  radiation  belts  at  any 
given  time.  1  herefore,  the  observations  are  a  constraint,  but  some  other  information  must  be  used  to 
“fill  in”  the  description  of  the  global  state.  In  physics-based  data  assimilation,  that  other  information 
is  the  “forecast”  state  provided  by  integrating  the  physical  equations  from  some  initial  condition  to  the 
current  time.  Then,  errors  in  the  physical  simulation  are  weighed  against  observational  errors  to 
produce  a  maximum  likelihood  estimate  of  the  actual  state  of  the  radiation  belts  the  “analysis”,  hi 
regions  where  there  is  no  data  or  where  the  data  is  of  low  quality,  the  forecast  prevails;  in  regions 
where  there  is  ample,  high  quality  data,  the  data  prevails;  elsewhere,  the  analysis  is  somewhere  in 
between  the  data  and  the  forecast.  In  other  words,  the  forecast  is  adjusted  as  little  as  possible  to  make 
it  consistent  with  the  observations. 

In  our  statistical  reanalysis,  the  role  of  the  forecast  is  replaced  by  the  empirical  statistical  model.  This 
“prior”  model  describes  the  likelihood  of  any  given  state  of  the  electron  belts;  thus,  we  choose  the  most 
likely  state  that  is  consistent  with  the  data  taken  in  a  window  around  each  6-hour  snapshot.  1  he  use 
of  a  statistical  prior  rather  than  physical  forecast  separates  the  “filling  in”  of  the  energy /spatial 
domain  from  specific  assumptions  about  physical  processes  (energization,  loss,  transport),  rhus,  the 
statistical  reanalysis  provides  a  unique  opportunity  to  test  physical  hypothesis.  1  he  statistical  model 
also  enables  Monte  Carlo  scenarios  that  are  realistic  and  portray  the  variety  of  possible  futures  a  space 
mission  might  encounter.  It  is  the  Monte  Carlo  scenario  capability  that  allows  our  model  to  provide 
error  bars  and  other  indications  of  uncertainty. 
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2  Statistical  prior  model 


Our  model  begins  with  an  a  priori  probability  model  of  the  statistical  variation  and  energy /spatial 
covariation  of  electron  fluxes.  This  “prior”  model  was  developed  from  historical  electron  observations. 
We  will  describe  the  mathematical  framework  for  the  prior  model  and  the  high  resolution  data  that 
was  used  to  develop  it. 

The  statistical  prior  model  does  not  include  an  explicit  solar  cycle  dependence.  Under  most 
circumstances,  such  a  dependence  is  unhelpful,  as  mission  schedules  slip,  and  one  would  not  wish  to 
have  to  redesign  the  radiation  protection  system  due  to  a  slip.  NASA-HDBK-4002  [1999,  page  22] 
states  “The  project  manager  knowing  the  mission  schedule  may  wish  to  assume  some  risk  in  order  to 
save  project  resources,  but  it  is  not  recommended  because  of  the  environmental  variation.” 

1  he  statistical  prior  model  includes  4  pieces:  a  spatial  map  of  the  distribution  of  fluxes,  an  error 
covariance  matrix  for  the  spatial  map,  a  spatial  covariance  matrix,  and  a  spat io temporal  covariance 
matrix.  In  this  context,  our  “spatial”  coordinates  are  energy  E ,  equatorial  pitch  angle  aeg,  and 
Mcllwain’s  Lm.  Throughout  this  work,  we  use  the  Olson-Pfitzer  Quiet  (OPQ)  magnetic  field  model 
[Olson  and  Pjitzer ,  1974]  added  to  the  IGRF  [IGRF,  2009],  as  this  combination  approximates  the 
“typical”  state  of  the  magnetic  field,  even  if  it  ignores  the  dynamics.  A  more  sophisticated  dynamic 
field  model  would  require  solar  wind  or  other  inputs  that  are  often  unknown  in  the  past  and  are 
certainly  unknown  for  future  missions. 

Our  energy  grid  spans  40  keV  to  10  MeV  in  12  logarithmically  spaced  points.  Our  pitch  angle  grid  has 
18  points  spanning  from  5  to  90  degrees  in  5  degree  steps.  Our  Lm  grid  spans  from  the  inner  edge  at 
~  1.13  to  11  in  27  points,  with  finer  resolution  (0.2)  in  the  inner  zone  gradually  increasing  to  coarser 
resolution  (1.0)  in  the  outer  zone.  For  many  computations,  we  wrork  on  a  reduced  grid  which  removes 
points  inside  the  loss  cone.  This  reduced  grid  has  418  of  the  possible  486  aeq  —  Lm  points.  The  total 
grid  size  is  5832  points,  which  reduces  to  Nx  =  5016  w  hen  the  loss  cone  points  are  removed. 

We  have  used  the  loss  cone  model  from  AE8  [Vette,  1991]  to  define  the  nominal  loss  cone.  Vette 
provides  an  empirical  formula  for  the  equatorial  pitch  angle  a/r  of  the  atmospheric  loss  cone: 

L  <  2.4 
2.4  <  L  <  3.0 
L  >  3.0 

At  Lm  =  1.1293  the  loss  cone  is  90  degrees  at  the  equator;  i.e.,  there  can  be  no  trapped  particles  belowr 
this  Lm  value.  Vette  used  an  IGRF-like  magnetic  field  model  to  define  Lm.  By  employing  the  AE8 
limits  to  Lm  computed  with  the  Olson-Pfitzer  model,  we  are  introducing  a  slight  inconsistency-one  we 
feel  confident  can  be  neglected  because  of  the  lowr  Lm  values  at  which  the  atmospheric  loss  cone  is 
substantial,  where  Olson-Pfitzer  adds  little  to  IGRF.  Evolution  of  the  IGRF  since  AE8  was  developed 
is  potentially  a  bigger  issue. 

The  map  of  the  distribution  of  fluxes  is  an  extension  of  the  traditional  map  of  the  mean  flux  on  a  grid. 
We  denote  by  x  differential  number  fluxes  of  electrons  on  our  model  grid.  We  adapt  the  flux  map 
concept  by  including  twro  parameters  of  the  statistical  distribution  of  flux  at  each  energy /spatial  grid 
point.  These  two  parameters  are 
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9i  =  logm50,  (2) 

02  =  log(mg5  -  m50),  (3) 

where  rn^ 0  and  mgs  are  the  50th  and  95th  percentiles  of  the  distribution  of  flux  at  each  grid  point.  The 
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inverse  is: 


m50  =  exp  (Oi),  (1) 

m95  =  exp(Oi)  +exp(02)-  (5) 

Throughout  this  paper,  we  use  log  to  represent  natural  logarithm  (log,,).  We  have  chosen  the  0 
parameterization  so  that  77195  >  W50  >  0  even  when  one  or  both  components  of  0  is  negative.  We  will 
explain  below  how  we  estimate  0  at  each  grid  point. 


We  will  need  to  track  the  estimation  errors  in  0  to  help  characterize  the  impact  of  measurement 
uncertainty  and  the  finite-duration  of  in  situ  observations  on  our  ability  to  specify  the  radiation 
environment  for  future  missions.  We  will  retain  the  error  in  that  estimation  process  in  an  error 
covariance  matrix  S2,  given  by: 


S  = 


S*  = 


-2  (kl) 


5(1) 

5(2) 

52(11)  52(12) 

52(21)  ^2(22) 

COV  (80k(Et, 

Bm,i  ),  80 1  ( Ej ,  <*eq,j  »  » 


(6) 

(7) 

(8) 


where  60  is  the  standard  error  estimated  for  0.  We  note  that  in  our  implementation.  5  is  not  a  square 
matrix,  and  details  of  its  calculation  will  be  given  below  in  section  2.2.6. 


We  will  need  the  spatial  covariance  matrix  to  ensure  that  the  radiation  belts  in  Monte  Carlo  scenarios 
and  in  the  reanalysis  are  appropriately  smooth.  L  he  spatial  covariance  matrix  is 

=  CO  v(zi,Zj),  (9) 

Zi  =  *-1(F(ar<;#<>))>  (10) 

where  F(x{:0^)  is  the  cumulative  marginal  distribution  for  x^ ,  and  4>  is  the  cumulative  distribution  of 
a  standard  Gaussian  (unit  normal).  I  he  operation  <t>  1  converts  Xi  from  a  variable  with 

the  desired  marginal  distribution  to  one  with  a  normal  distribution.  We  have  shown  in  [O'Brien,  2005] 
that  the  Weibull  distribution  is  an  excellent  representation  of  the  marginal  distribution  of  electron 
fluxes  at  fixed  grid  points. 


Figure  1  compares  the  observed  distribution  in  a  selected  bin  to  the  Weibull  distribution.  In  the  top 
panel,  a  quantile-quantile  plot  shows  that,  with  the  exception  of  the  lowest  fluxes  in  the  distribution, 
the  Weibull  is  a  good  representation  of  the  distribution.  1  he  bottom  panel  shows  a  similar  result  in 
the  form  of  the  cumulative  distribution  (P<,  or  F)  and  exceedance  probability  (P>,  or  1  —  F). 


I  he  transformation  from  x  to  z  allows  us  to  utilize  the  wide  array  of  tools  for  multivariate  Gaussians 
to  model  the  spatial  and  temporal  smoothness  of  the  radiation  belts.  Formally,  we  are  using  a 
Gaussian  copula  [see,  e.g.,  Nelsen ,  1999],  which  describes  the  relationship  between  Xi  and  Xj  by  way  of 
transformed  variables  (ui,uj),  where  U{  =  F(xi]0^),  which  makes  u  equivalent  to  the  fraction  of  the 
data  that  falls  below  x.  Both  ut  and  Uj  are  uniformly  distributed  from  0  to  1,  but  they  retain  the 
correlation  information  for  xt,xj.  By  transforming  to  uniform  variables  in  the  [0,  1]  domain,  we  can 
model  the  correlation  structure  independently  of  the  marginal  distributions.  Figure  2  shows  that  the 
scatter  plot  of  transformed  daily  average  flux  in  two  bins  clusters  in  the  regions  of  high  probability 
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Quantiles  of  Weibull,  #/cm2/s/sr/keV 
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Figure  1:  An  example  of  the  Weibull  distribution  fit  to  observed  fluxes  in  a  selected  energy  channel  in  a 
bin  of  aeq  and  Lm.  The  distribution  fits  well  everywhere  except  the  very  lowest  2-3%  of  the  distribution, 
which  is  of  little  consequence  for  spacecraft  effects.  The  top  panel  is  a  quantile-quantile  plot,  and  the 
bottom  panel  show's  the  tails  of  the  cumulative  distribution  and  exceedance  probability. 
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Gaussian  Copula  forr=0.71,  c=0:0.5:4 
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u,  1368  keV,  n  =40°  L  =4.7 

’  eq  m 


Figure  2:  A  pairwise  scatter  plot  of  the  transformed  fluxes,  u,  for  two  bins,  superimposed  on  the 
probability  density  of  the  Gaussian  copula  with  correlation  parameter  0.71  taken  from  the  data.  Contour 
lines  are  provided  at  copula  densities  of  0  to  4  in  steps  of  0.5.  As  expected,  the  points  cluster  where  the 
density  is  high  (red). 
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indicated  by  the  copula.  The  copula  has  a  single  parameter  for  each  pair  of  grid  points,  which  is  given 
by  rij  =  Ejj,  our  spatial  correlation  coefficient  for  the  z's  given  in  (9). 

The  Weibull  distribution  is  given  by: 


F(x)  = 

1  -exp[-(>/x0)7], 

(11) 

7  = 

'"e  [SfrSf  ] 

(12) 

log  ^ 

G  77150 

Xq 

m50 

(13) 

[-  log(l  —  0.50)]1/7 

For  reference,  the  average  value  of  a  Weibull  distribution  is  given  by  XoF(l  4-  I/7),  where  T(c)  is  the 
complete  Gamma  function  with  argument  c. 

A  useful  alternative  to  the  Weibull  distribution  (one  that  seems  appropriate  for  the  inner  zone  proton 
belt),  is  the  log-normal,  given  by: 


F(x) 

-  *(^)-)' 

(14) 

=  logmso, 

(15) 

CT 

log  mgS  -  fl 

(16) 

For  reference,  the  mean  of  a  log  normal  distribution  is  given  by  exp  (/z  -f  <t2/2) 

Finally,  in  order  to  have  appropriate  smoothness  in  time,  our  model  uses  a  spatiotemporal  covariance 
matrix  R ,  which  is  given  by 


Rij  =  co v(zi(t),Zi(t  *4*  St))  ,  (17) 

where  St  is  1  day.  R  gives  the  correlation  at  one  day  lag,  which  means  the  model  will  not  effectively 
describe  the  variation  in  the  radiation  belts  on  timescales  shorter  than  a  day.  This  is  largely  a 
limitation  of  the  observational  sampling  in  and  «e<7,  which  are  tied  to  the  orbital  period  of 
satellites  and  the  spin  period  of  the  Earth. 

The  matrices  E  and  R  will  be  processed  into  Q,  G,  and  C,  which  are,  respectively,  the  principal 
components  matrix  of  E,  the  6-hour  persistence  operator  for  the  principal  components,  and  the  noise 
conditioning  matrix  for  the  principal  components.  The  details  of  these  latter  matrices  will  be  given  in 
section  2.3.  Together,  these  data  tables,  along  with  6  and  5  on  the  model  grid,  describe,  to  first  order, 
the  statistical  properties  of  the  electron  belts,  along  with  an  indication  of  the  model  uncertainty. 


2.1  Data  used  to  construct  the  statistical  prior  model 

To  build  our  statistical  prior  model,  we  have  used  data  from  CURES  MP]A  [  Vampola  et  a/.,  1992], 
CRRES  HEEF  [Dichter  et  a/.,  1993],  S3-3  MES  [Vampola  and  Adams ,  1988],  SCATHA  SC-3  [Reagan  et 
a/.,  1981],  and  Polar  CEPPAI)  [Blake  et  at,  1993].  All  of  these  instruments,  to  some  degree,  have 
calibration  errors  and  backgrounds  that  ought  to  be  removed.  At  time  of  writing,  that  process  is  not 
complete.  Instead,  we  have  opted  to  combine  all  of  the  data  with  nominal  calibrations,  and  assume 
that  the  remaining  errors  and  backgrounds  will  be  different  enough  among  the  platforms  that  they  will 
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Table  1:  Data  used  to  build  the  statistical  prior  model 
Vehicle  interval  Instruments  Energy  Coverage  L 


S3-3 

1976-1979 

MES 

0.01-1.6  MeV 

all 

SCATHA 

1979-1990 

SC-3 

0.05-4.6  MeV 

5-8 

CRRES 

1991-1992 

MEA/HEEF 

0.1-7  MeV 

2.5-7 

Polar 

1996-2008 

CEPPAD 

0.02-10  MeV 

all 

be  smoothed  out  when  we  combine  the  data  from  multiple  sensors.  This  assumption  has  largely  been 
proven  true.  For  example,  Figure  10  shows  a  spectral  lit  of  CURES  and  SCATHA  data  in  one 
aeq  “  Lm  bin.  At  high  energies,  the  SCArl  HA  background  causes  its  spectrum  to  curve  upward, 
unphysically.  However,  the  nearest- neighbors  interpolation  and  smoothing  algorithm  smoothes  out  this 
feature  (see  section  2.2  for  details). 

Table  1  provides  an  overview  of  the  time,  energy,  and  L  range  of  the  data  used  to  build  the  prior 
model.  These  data  were  selected  because  they  are  energy-  and  angle- resolved  measurements,  so  that  no 
sophisticated  inversions  are  required  to  convert  them  to  the  unidirectional  differential  fluxes  need  for  x 
in  the  model. 

We  have  filtered  out  some  data  which  is  known  to  be  suspect:  CRRES/MEA  data  for  Lm  <  3,  due  to 
possible  proton  and  bremsstrahlung  background:  the  ~  400  keV  channel  on  Polar/CEPPAD  was 
removed  as  it  is  often  much  higher  than  neighboring  points  in  the  merged  IES/HlSTe  spectrum; 
S3-3/MES  channels  1-4  were  removed  at  the  suggestion  of  Vampola’s  data  release  notes.  Known  issues 
that  were  not  addressed  include:  cosmic  ray  backgrounds  in  SCATHA/SC-3,  dead-time  in 
Polar/CEPPAD/HISTe,  and  ion  contamination  in  Polar/CEPPAD/IES.  1  hese  issues  will  be  addressed 
in  later  models  (i.e.,  AE9). 


2.2  Data  processing  to  develop  the  prior  model 

Figure  3  provides  an  overview  of  the  data  processing  steps  needed  to  obtain  the  model  data  tables,  0 , 
g,  £,  and  Q. 

The  data  sets  begin  as  daily  files  of  unidirectional  differential  flux  (#/cm2/sr/s/keV)  as  a  function  of 
time.  The  data  are  provided  either  as  energy  spectra  versus  time,  x(Ek,&ty  Tm>t),  or  energy-angle 
spectra  versus  time,  x(Ek,oti ,  Lm,t)-  ln  either  case,  the  pitch  angle  a  is  a  local  pitch  angle,  and  Lm  is 
given  for  locally  mirroring  particles  only.  That  is,  the  angular  dependence  of  Lm  is  ignored,  which 
introduces  a  minor  inconsistency,  but  one  no  worse  than  the  assumption  that  Lm  and  aeq  are  invariant 
over  a  particle  drift  ( aeq  is  a  bounce  invariant  only).  Because  at  this  time  we  have  no  error  estimates 
for  our  data,  we  assume  a  fixed  10%  error  for  all  fluxes:  Slogx  =  0.1. 


2.2.1  Binning 

The  first  step  is  to  assign  the  time-series  data  to  aeq  —  Lm  bins  so  that  we  can  compute  statistics  over 
multiple  time  samples  at  the  same  phase-space  location.  Because  we  do  not  know  what  the  best 
aeq  ~  Lm  grid  is  to  begin  with,  we  use  a  dynamic  binning  approach  to  organize  each  data  set.  We  bin 
each  data  set  separately.  The  dynamic  binning  recursively  bisects  the  aeq  and  Lm  domain  to  form 
contiguous  bins  that  meet  the  following  criteria: 
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Multiple  Sensors 


Figure  3:  Flow  chart  depicting  data  processing  steps  that  lead  to  the  prior  model  data  tables. 
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•  1  he  bin  spans  at  least  1  degree  in  aeq  and  0.05  in  Lm, 

•  The  first  and  last  sample  in  the  bin  are  separated  by  at  least  182.5  days  (1/2  year), 

•  The  bin  includes  at  least  100  unique  days  of  data. 


Once  the  bisection  is  complete,  any  bin  that  spans  more  than  5  degrees  in  aeq  or  0.25  in  Lm  is 
excluded  from  further  processing.  I  his  dynamic  binning  scheme  allows  us  to  optimize  the  aeq  —  Lm 
resolution  of  the  flux  maps  relative  to  the  available  data  from  each  mission.  We  can  chose  the 
resolution  of  the  model  grid  after  we  have  already  binned  the  in  situ  data. 

We  note  that  the  energy  dimension  of  the  bins  is  slightly  different  from  the  aeq  and  Lm  dimensions: 
each  sensor  provides  fluxes  on  a  unique,  but  pre-determined,  energy  grid  (the  energy’  channels  /J*.)- 
Therefore,  each  sensor  on  each  spacecraft  has  its  own  idiosyncratic  set  of  bins  in  E  as  well  as  in  c\eq, 
and  Lm. 

When  the  binning  is  complete,  each  time  (and,  possibly,  angle)  sample  has  been  assigned  to  a 
cteq  —  Lm  bin.  Figure  4  shows  the  number  of  days  per  bin,  in  selected  energy  ranges,  combined  across 
all  the  missions  in  table  1.  CR11ES  provides  a  minimum  of  hundreds  of  days  over  nearly  the  entire 
domain  (see  also  Figure  5).  SCATHA  provides  a  broad  band  of  long-term  coverage  near  GEO  at  all 
pitch  angles.  Polar  provides  long-term  coverage  at  lower  pitch  angles. 


2.2.2  Daily  averaging  within  bins 

Next,  we  compute  daily  averages  and  standard  deviations  for  each  bin,  in  a  two-step  process.  We  do 
this  separately  for  each  sensor  on  each  vehicle.  First,  we  compute  the  average  for  each  pass  through 
the  bin  during  the  day,  then  we  compute  the  average  of  those  passes  to  give  a  daily  average  in  the 
energy /spatial  bin.  A  pass  is  defined  as  a  set  of  time  samples  in  the  bin  with  no  gaps  longer  than  15 
minutes.  In  this  way,  a  long  pass  and  a  short  pass  through  the  same  bin  would  have  the  same  weight  in 
the  daily  average.  1  his  minimizes  the  impact  of  serial  correlation  in  the  time  series  data. 

The  averages  are  computed  with  weighting  based  on  the  estimated  relative  error  (rr)  in  the  observed 
fluxes,  as  follows: 


<  x  >  = 


9 

<  X*  > 


5>r2’ 

l>r2’ 


<  x 


—  1. 


(18) 

(19) 

(20) 


In  the  first  averaging  step  (pass-averaging),  the  weights,  rrt  are  simply  the  d logx  from  the  original  time 
series.  In  the  second  step  (daily- averaging),  they  are  given  by  the  results  of  (20)  in  the  previous 
(pass- averaging)  step.  If,  after  either  step  of  the  averaging  process,  (t<x>  is  less  than  0.1,  it  is  replaced 
by  0.1.  Similarly,  if  it  is  greater  than  log  10,  it  is  replaced  by  log  10.  1  liese  ad  hoc  limits  ensure  that 
the  data  analysis  process  does  not  idiosyncratically  introduce  unrealistically  small  or  unrealistically 
large  error  bars. 
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Figure  4:  TEM-2  data  coverage  combined  over  all  missions  in  selected  energy  ranges. 
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2.2.3  Computing  0  and  cov 60  in  each  bin 


With  the  data  binned  and  daily  averaged,  we  compute  0  for  each  bin  for  each  sensor  on  each  vehicle. 
Taking  all  the  data  for  a  particular  energy /spatial  bin,  we  perforin  a  bootstrap  resampling  operation  to 
estimate  the  error  covariance  co v60  in  the  bin.  The  bootstrap  amounts  to  randomly  resampling  the 
array  of  daily  averages,  and  then  perturbing  each  selected  daily  average  by  a  random  multiplicative 
factor  exp[(T<x>c],  where  e  is  a  random  number  chosen  from  a  standard  normal  distribution.  1  his 
resampling  procedure  gives  us  an  initial  local  estimate  of  cov 60,  including  the  covariance  of  errors 
between  0\  and  02. 

We  adjust  this  initial  estimate  in  two  ways.  First,  we  inflate  the  error  covariance  matrix,  if  necessary, 
to  account  for  serial  correlation  among  the  daily  averages.  Let  p\  be  the  linear  correlation  between 
adjacent  daily  averages  in  the  bin.  The  error  covariance  should  be  multiplied  by  (1  +  pi)/(l  —  p\) 

[ Zwiers  and  von  Storch ,  1995],  as  a  first-order  estimate  of  the  effect  of  serial  correlation.  Then,  this 
second  estimate  of  the  error  covariance  is  rescaled  as  needed  to  ensure  that  the  variance  of  60\  and  602 
are  both  at  least  0.12  (i.e.,  10%  relative  error).  Conversely,  if  the  variance  of  either  60k  is  larger  than  1 
(100%  relative  error),  the  associated  value  of  0  is  omitted  from  further  analysis  because  the  error  is  too 
large.  The  variance  rescaling  conserves  the  linear  correlation  coefficient  implied  by  the  diagonal  term  in 
the  covariance  matrix.  We  denote  the  adjusted  covariance  matrix  in  bin  i  with  s^: 

^ki  ~  ^ov  (60fc(/yj,  Lm,t),  60/(/>i,  oieg  j,  Lm>j))  (^1) 


Figure  5  shows  0,  errors,  error  correlation  ^r(60i,602)  =  s i 2* /  \  s\ V ,s 22*  J  »  and  data  coverage  for  one 

CRRES  MEA  channel.  The  lack  of  data  at  high  Lm  arises  from  the  CRRES  orbit,  and  the  lack  of  data 
at  low  Lm  arises  from  the  removal  of  Lm  <  3  data,  as  indicated  in  section  2.1.  Also,  because  the  orbit 
was  inclined,  coverage  at  high  Lm  tended  to  be  ofF  the  equator,  producing  a  lack  of  data  at  high  aeq  at 
high  Lrn.  Estimated  uncertainty  in  0  is  largest  at  low  Lm,  where  variability  is  most  dramatic,  but  the 
correlation  between  errors  in  0\  and  02  is  small. 


2.2.4  Interpolating  to  the  model  grid 

Next,  we  combine  all  of  the  0  from  all  of  the  bins  from  all  of  the  sensors  into  one  global  table  of 
0(Ei,aeq>i,  Lm,i)  on  the  model  grid. 

First,  we  fit  a  “base”  model  of  0  by  the  method  of  least  squares: 

o\0)(E,aeq,Lm))  =  n4  logsin  max  {Qefl -ojc(£-m),  0.01] 

+  log  E  +  Ci(Lm  -  4.5)2.  (22) 

This  model  amounts  to  a  power  law  in  energy,  a  roughly  sinn  dependence  on  equatorial  pitch  angle,  a 
single  belt  centered  at  Lm  =  4.5,  and  a  loss  cone  given  by  (1).  The  coefficients  are  given  in  Table  2. 

1b  interpolate  to  the  model  grid,  we  have  taken  the  average  of  120  nearest  neighbors.  The  nearest 
neighbors  are  defined  in  terms  of  the  relative  distance  from  the  grid  point  to  the  edge  of  the  domain. 
Equatorial  pitch  angle  and  Lm  are  linearly  transformed  onto  the  [0,  1]  domain,  and  energy  is 
transformed  logarithmically.  Distance  in  the  energy  dimension  is  further  scaled  down  by  a  factor  of  2 
to  ensure  smooth  spectra. 
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CRRES  MEA,  876  keV 


Figure  5:  Median  and  95th  percentile  flux  maps  for  one  CRRES  MEA  channel,  along  with  errors  on  0 
and  coverage.  A  black  “X”  indicates  that  an  aeq  —  Lm  bin  did  not  meet  the  bin  selection  criteria  defined 
in  section  2.2.1,  usually  because  of  too  few  days  of  data. 


Table  2:  Coefficients  of  0^ 

%  CLi  b-i  flj  Cj 

1  206  2J38  0187  -0.0684 

2  21.3  2.41  0174  -0.0728 
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We  will  need  the  inverse  of  at  each  grid  point,  but  we  will  organize  this  inverse  as  a  large  matrix 
that  spans  the  energy/spatial  grid  twice  on  each  side  rather  than  as  a  set  of  small  2  x  2  matrices: 

<!?’=<>«  ([s‘T')w.  <*» 

where  Sij  is  the  Kronecker  (discrete)  delta  function. 

At  each  grid  point  (?'),  we  find  the  120  nearest  neighbors  (j),  and  compute  a  weighted  average  based  on 
weights  provided  in  s  via 


^  kj  — ®k ( i  ®eq,j 5  @  ^  (  Ey i  &eqji  ^mj)) 


£<n>r  £<i2>r 

^(21)  f  ^(22)  J*  - 

■  c(H)x+c(l2)x " 

£(2i)y  _j_  ^(22)y  ’ 


'  A(»  ' 
A<2> 


'I'-'T 


(24) 

(25) 

(26) 

(27) 


0k(Ei,cte<l'i,  Lm,i)  -  A<fc)  +e{°)(Ei,  aeq,i,  Lm,i). 


(28) 


After  completing  this  procedure,  we  have  an  estimate  of  0  at  each  grid  point  (as  opposed  to  our  earlier 
0  estimates,  which  were  organized  in  non-uniform  bins  unique  to  each  sensor).  We  call  this  grid  of  0  a 
flux  map  because  it  can  be  used  to  reconstruct  the  median  and  95th  percentile  fluxes  on  the  model  grid. 


2.2.5  Flux  map  diagnostics 


In  this  section,  we  provide  several  figures  that  characterize  the  flux  maps  ( 0 ).  Figure  6  provides 
aeq  —  Ijm  flux  maps  of  the  median  and  95th  percentile  fluxes  derived  from  6  according  to  (4)  and  (5)  at 
several  fiducial  energies.  The  combined  flux  map  is  typically  not  as  smooth  as  the  flux  maps  from 
individual  sensors  (e.g.,  compare  to  Figure  5).  1  herefore,  in  the  subsequent  development  of  of  AF9  and 
AR9,  we  expect  to  revise  the  nearest -neighbors  algorithm. 

figure  7  shows  how  TEM-2  compares  to  AE8  and  CRRESELE.  We  have  plotted  the  median  from 
1  EM-2  to  enable  comparison  to  AE8,  which  is  actually  closer  to  a  median  flux  model  than  it  is  to  an 
average  flux  model  [  Vette,  1991]  (the  AE8  flux  map  results  from  a  least-squares  fitting  to  observed  log 
fluxes,  which,  for  the  log-normal  distribution  used  by  Vette ,  results  in  a  fit  to  the  median  flux). 
CRRESELE  does  not  provide  a  median  model,  but  it  does  provide  an  average  [Brautigam  and  Bell , 
1995],  which  allows  a  rough  comparison  to  AE8  and  TEM-2  medians.  What  we  see  is  that,  as 
expected,  the  models  do  not  agree  quantitatively.  Surprisingly,  CRRESELE  places  the  outer  belt  at 
systematically  lower  radial  distance  than  does  AE8  or  TEM-2.  That  is  especially  surprising  considering 
that  TEM-2  relies  heavily  on  CRRES  data.  At  1.1  MeV,  all  of  the  models  are  in  considerable 
disagreement  as  to  the  location  of  the  major  topological  features  (inner  belt,  slot,  outer  belt). 

Figure  8  shows  the  same  three  models  as  in  Figure  7  but  along  a  field  line  through  the  outer  zone. 
Again,  disagreement  is  profound.  Generally,  it  implies  that  TEM-2  pitch  angle  distributions  are  more 
isotropic  (flatter)  than  those  implied  in  AE8.  CRRESELE  appears  to  be  the  most  isotropic  (flattest). 
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Median  Flux  95%  Flux 


Figure  6:  TEM-2  maps  of  median  and  95th  percentile  fluxes. 


Figure  7:  Comparison  of  TEM-2  to  AE8  and  CRRESELE  along  an  equatorial  radial  profile. 
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Latitude  Profile  in  TEM-2  Median  Along 
Field  Line  Through  4. 5, 0,0  MAG  at  1 -Jan-2000 
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Figure  8:  Comparison  of  TEM-2  to  AE8  and  CRRESELE  along  a  field  line  through  (4. 5, 0,0)  Re  hi  tht 
outer  zone. 
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Latitude  (Geodetic) 

Figure  9:  Comparison  of  TEM-2  to  AES  and  CRRESELE  along  a  LEO  trajectory  through  the  South 
Atlantic  Anomaly. 


Figure  9  shows  the  same  three  models  as  in  Figure  7  but  along  a  low  altitude  longitude  line  through 
the  South  Atlantic  Anomaly.  Yet  again,  disagreement  is  profound.  Even  though  TEM-2  and  AE8 
parameterized  their  loss  cones  the  same  way  in  terms  of  Lm,  it  appears  that  the  differences  between 
IGRF  and  OPQ  cannot  be  neglected  when  defining  the  maximum  ratio  of  the  mirror  to  equatorial  field 
strength.  There  is  also  the  possibility  that  the  naive  treatment  of  the  angular  response  of  the  sensor 
data  has  led  to  an  artificially  full  loss  cone.  Specifically,  a  wide-angle  sensor  can  still  observe  particles 
in  its  aperture  when  looking  directly  along  the  field  line,  giving  the  incorrect  impression  of  particles  in 
the  loss  cone. 


2.2.6 


Computing  the  errors  in  the  interpolated  flux  map  0 


The  error  covariance  S2  for  the  0  flux  map  is  computed  in  a  bootstrap  fashion,  where  we  bootstrap  the 

6  flux  map  calculation  over  random  resamples  of  the  set  of  missions  from  which  the  0  flux  map  was 
derived.  Our  assumption  is  that  each  mission  may  have  its  own  residual  systematic  errors  in  the  flux 
map  it  would  produce  in  isolation.  Those  errors  arise  from  calibration  differences  between  sensors  on 
different  missions  and  from  differences  in  the  natural  environment  experienced  by  each  mission.  By 
bootstrapping  over  missions,  we  include  in  5  a  spatially-cor related  measurement  error  in  addition  to 
the  local  measurement  error  in  (21). 
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In  each  bootstrap,  we  also  randomly  perturb  each  0{l\  assuming  the  errors  have  zero  mean  with 
covariance  given  by  This  perturbation  propagates  the  local  error  (finite  sample  size,  measurement 
residual)  in  through  the  bootstrap,  into  5. 


We  compute  S_  directly  from  N^s 


=  50  bootstrap  solutions  for  0: 


e(fc) 


(k,m) 


(29) 

(30) 


where  m  identifies  the  mth  bootstrap  sample.  With  this  definition,  S_ S_ 1  is  the  maximum  likelihood 
estimate  of  the  covariance  matrix  for  6  given  the  bootstrap  samples. 


Figure  10  shows  the  spectrum  of  median  and  95th  percentile  fluxes  in  one  aeq  —  Lm  bin  for  three 
spacecraft  and  for  the  nearest-neighbors  fit,  with  error  bars.  For  the  most  part,  as  expected,  the  data 
lies  near  the  fit.  The  spread  of  the  data  is  not  large  in  this  bin,  and  the  various  measurements  all  agree 
fairly  well.  1  his  is  not  always  the  case,  but  documenting  the  variety  of  types  of  disagreement  between 
the  measurements  is  premature,  as  the  data  have  not  yet  been  intercalibrated. 


2.2.7  Computing  spatial  and  spatiotemporal  covariances 


The  calculation  of  spatial  (£)  anti  spatiotemporal  (R)  covariances  poses  a  significant  computational 
problem.  We  have  elected  to  calculate  a  limited  set  of  covariances  between  daily-averaged  fluxes  in  bins 
from  different  sensors,  and  then  to  interpolate  that  set  of  covariances  onto  the  model  grid  by  taking  the 
average  of  nearest  neighbors.  We  have  chosen  to  compute  3  x  105  spatial  covariances  and  twice  that 
many  spatiotemporal  covariances.  We  require  at  least  100  common  days  (offset  by  1  when  computing 
the  spatiotemporal  covariance)  between  any  two  bins  in  order  to  compute  a  valid  covariance. 


According  to  (9)  and  (17),  we  must  compute  the  covariances  on  Gaussianized  variables  £,  which  are 
transformed  from  the  observed  fluxes.  A  member  of  the  set  of  observed  daily  average  fluxes  in  a  bin  is 
denoted  x*.  We  denote  the  position  xt  would  have  in  a  sorted  set  from  least  to  greatest  as  ji ,  with 
positions  indexed  from  1  to  N.  We  can  then  empirically  convert  from  X{  to  Z{  as  follows: 


Zi  =  <*>-' 


N  +  1 


(31) 


This  empirical  formula  allows  us  to  convert  observed  fluxes  to  Gaussian  variables  without  explicitly 
assuming  a  marginal  distribution  (e.g.,  Weibull).  1  he  covariance  is  then  computed  between  the  z  from 
one  bin  and  the  z  from  another  bin. 


To  interpolate  to  the  model  grid,  we  have  taken  the  average  of  10  nearest  neighbors,  using  the  same 
transformed  coordinates  we  used  for  0 ,  but  without  the  energy  dimension  scaling  factor,  and  without 
any  error  weighting.  We  note  that  for  covariance  calculations  the  nearest  neighbors  are  in  a  6 
dimensional  space.  Although  there  is  considerable  error  in  the  computed  spatial  and  spatiotemporal 
covariances,  we  neglect  it  in  further  calculations  because  we  only  need  a  first-order  model  of  spatial 
and  spatiotemporal  covariance. 
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Median  Flux,  aeq=75,Lm=5.5 


c/) 


Figure  10:  Median  (top)  and  95th  percentile  (bottom)  fluxes.  Data  from  in  situ  observations  in  the  bin 
are  plotted  in  color.  The  nearest  neighbors  fit  is  indicated  in  black.  Also  shown  is  the  “base  0”  given 
by  22.  Error  bars  on  the  points  and  the  fit  indicate  one  standard  deviation. 
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2.2.8  Principal  components 


I  he  computed  spatial  covariance  mat  rix  E  is  singular  and  typically  not  positive  definite,  as  is  common 
for  large  covariance  matrices,  especially  those  constructed  from  the  kind  of  sparse  spatial  coverage  we 
have  in  the  radiation  belts.  Subsequent  operations  in  the  development  and  use  of  TEM-2  imply  the 
inverse  of  E.  To  avoid  explicitly  computing  E_1  we  employ  a  linear  change  (and  reduction)  of  variables 
from  z  to  q  that  will  represent  most  of  the  variation  in  the  z’s  while  having  a  non-singular,  positive 
definite  covariance  matrix. 

The  transform  will  have  the  following  form: 


z  =  (32) 

To  construct  Q ,  we  first  compute  the  eigen-values  and  eigen- vectors  of  E: 

(zzr)  =  £  =  «  Qgr.  (33) 

The  columns  of  V  are  orthonorinal,  and  they  represent  the  principal  components  of  variation  of  the 
transformed  fluxes  z.  D  is  a  diagonal  eigenvalue  matrix  with  nonincreasing  values  dt  along  the 
diagonal.  Wilks  [1995,  pp.  380-383]  reviews  a  number  of  methods  for  selecting  how  many  principal 
components  to  retain.  All  of  those  heuristics  selected  principal  components  of  E  that  are  clearly  noise. 
Therefore,  in  a  somewhat  ad-hoc  manner,  we  select  all  the  principal  components  that  explain  more 
that  1%  of  the  variance  over  all  gridpoints  (i.e.,  Dtl  >  Nx/ 100).  The  other  components  are,  we 
presume,  noise. 

Thus  we  construct  V  and  D  by  retaining  only  the  first  J\rq  =  8  columns  of  V  and  the  first  Nq  rows  and 
columns  of  I) .  Our  initial  estimate  of  Q  is  then 


g  =  £|>l/2-  (34) 

*  a  T 

However,  because  we  have  left  out  some  eigenvalues,  we  know  that  QQ  only  approximates  E-  We  will 
need  a  refined  estimate  of  Q  that  preserves  unit  variance  for  all  the  z’s.  'That  is,  we  need  to  scale  Q  so 
that  the  diagonal  elements  of  QQ1  are  one. 

The  correction  matrix  follows  from: 

Aij  =  Stj(Qg  )u.  (35) 


We  can  now  define  Q  and  its  psendoinverse  as  Qf  as: 

(36) 

(37) 


Q<  =  d-i,2vta,,\ 


Thus  we  have  avoided  the  problem  that  E  is  singular  by  eliminating  the  redundant  variables.  The  role 
of  E'1  is  now  played  implicitly  by  E*  =  ( Q i)1  Qi . 


Because  the  transform  from  q  to  z  is  linear  with  no  additive  offset,  the  q's  are  Gaussians  with  mean 
zero.  Because  the  q' s  are  the  coefficients  of  the  eigenvectors  of  E,  the  covariance  matrix  for  the  q  s  is 
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simply  /;  that  is,  the  q's  are  linearly  independent  unit  Gaussians.  The  q' s  represent  the  amplitudes  of 
the  principal  independent,  normalized  components  of  the  energy /spatial  structure  of  the  electron  belts. 
They  are  loosely  analogous  to  the  principal  Fourier  components  of  a  simple  harmonic  oscillator  after  a 
low-pass  filter. 


2.3  Engineering  application:  Monte  Carlo  mission  scenarios 

If  we  can  generate  a  realistic  surrogate  time  series  of  qt,  then  we  can  transform  back  into  fluxes 
q  — *  z  — *  x  to  obtain  a  realistic  surrogate  time  series  of  global  fluxes.  Such  a  surrogate  time  series 
retains  the  statistical  distribution  and  spatial  covariances  of  the  fluxes,  and  can  be  used  as  a  Monte 
Carlo  mission  scenario  for  engineering. 

In  order  to  generate  a  realistic  series  qt ,  we  must  use  R  and  Q  to  construct  a  first-order  autoregressive 
(AR1)  process  for  the  q's  that  will  yield  the  right  spatiotemporal  covariance  of  the  fluxes,  as  well.  An 
AR1  process  for  q  can  be  written  as: 

qt+&t  —  G.qt  +  (38) 

where  fjt  is  uncorrelated  standard  Gaussian  white  noise,  and  we  initialize  qo  to  uncorrelated  standard 
Gaussian  white  noise. 

With  some  matrix  manipulations,  we  can  obtain: 


( ;  = 

rr>-l  St fT 

s' if  (s') 

(39) 

£  = 

(40) 

The  matrix  powers  are  taken  via  eigenvalue  decompositions.  During  this  operation,  the  eigenvalues  of 
G  are  limited  to  lifetimes  of  3  years,  and  imaginary  eigenvalues  of  C  are  set  to  zero.  These  numerical 
corrections  prevent  unrealistically  long-lived  or  growing  eigenmodes  and  imaginary  values.  As  it  turns 
out,  in  the  final  build  of  TEM-2,  neither  of  these  corrections  was  needed. 

Using  equation  (38),  we  can  generate  an  arbitrary  number  of  realistic  scenarios  for  qt  and,  in  turn,  xt. 
That  is,  we  can  generate  an  arbitrary  number  of  realistic  Monte  Carlo  scenarios  for  the  time  series  of 
the  global  flux  environment  over  the  mission.  It  is  then  a  relatively  simple  matter  to  compute  from 
each  time-dependent  global  flux  scenario  the  time  series  of  flux  at  a  notional  spacecraft  as  it  moves 
along  its  orbit  over  a  mission.  For  each  scenario,  we  can  compute  desired  engineering  quantities,  such 
as  total  mission  dose  versus  depth,  or  worst  case  T-hour  fluence. 

In  order  to  capture  the  measurement  error,  i.e.,  our  uncertainty  in  the  underlying  flux  map  6 ,  we 
perturb  the  0  to  a  unique  0  map  for  each  scenario,  according  to: 

-^m,i)  =  ^fc(-^i,  &eq,i > 

+E5S)ei-  (4l) 

3 

where  each  tj  is  a  random  number  chosen  independently  from  a  standard  Gaussian  distribution,  and  S 
is  given  by  (6). 

Thus,  even  in  long  missions,  for  which  dynamic  variability  may  average  out,  each  scenario  will  give  a 
slightly  different  effect  (e.g.,  total  dose)  consistent  with  our  uncertainty  in  0. 
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>2  MeV  Electrons  at  136°  W  Longitude,  TEM-2  and  GOES-10 


Day  of  March  2006 

Figure  11:  Comparison  of  TEM-2  scenarios  and  reanalysis  to  AE8  and  GOES-10  observations  for  a  real 
interval.  GOES-10  flux  is  multiplied  by  3 tt  to  approximate  conversion  from  unidirectional  to  omnidirec¬ 
tional  flux,  assuming  that  the  flux  follows  a  sin5/4  a  distribution  in  local  pitch  angle.  CRRESELE,  as 
plotted,  is  driven  by  “Apl5” ,  the  15-day  trailing  average  of  the  Ap  index  [Brautigam  and  Bell ,  1995]. 
The  percentiles  shown  are  taken  from  200  Monte  Carlo  scenarios,  of  which  only  20  are  shown  (gray). 


2.3.1  GEO  example 

Figure  11  shows  a  short  interval  of  a  comparison  between  TEM-2  (including  several  Monte  Carlo 
scenarios)  and  observed  GOES  >  2  MeV  electron  flux.  As  expected,  quantitative  agreement  is  poor  in 
fact,  GOES-10  spends  a  significant  amount  of  time  measuring  fluxes  well  below  the  5th  percentile 
computed  from  200  scenarios.  However,  the  different  TEM-2  scenarios  show  some  of  t  he  diurnal 
variation  evident  in  the  GOES  data,  and  they  clearly  show  dynamic  behavior,  which  is  absent  from 
AE8.  1  lie  very  short  timescale  jagged  oscillations  on  the  timescale  of  an  hour  or  so  are  likely  numerical 
artifacts  of  the  calculation  of  the  magnetic  coordinates. 

What  the  TEM-2  prior  model  is  really  built  to  do  is  represent  the  long-term  distribution  of  electron  flux 
along  any  orbit.  Figure  12  shows  a  comparison  of  statistical  distributions  from  TEM-2  scenarios  versus 
GOES-10  for  >  2  MeV  electron  flux.  The  format  is  the  same  as  Figure  1.  The  top  panel  shows  that 
the  distribution  produced  by  TEM-2  can  be  off  by  as  much  as  a  factor  of  10  from  that  of  GOES- 10, 
and  it  does  not,  generally  have  the  same  shape.  1  he  only  part  of  the  distribution  that  TEM-2 
represents  well  is  the  highest  5%.  Preliminary  results  with  the  AE9  beta  reveal  that  the  shape  problem 
will  be  resolved,  but  there  will  be  some  residual  ofTset  (factor  of  ~  3).  Therefore,  we  can  att  ribute  the 
shape  discrepancy  to  our  choice  of  magnetic  coordinates  and  the  absence  of  GEO  data  in  TEM-2. 

Extending  the  simulation  shown  in  Figure  11  to  an  ent  ire  10-year  interval,  Figure  13  provides  spectra 
of  means  and  24-hour  (exponentially-weighted)  worst  cases  for  various  models,  including  TEM-2. 

These  worst-case  spectra  are  an  approximation  to  what  one  would  need  for  assessing  an  internal 
charging  hazard.  In  practice  it  would  be  better  to  run  an  actual  time-dependent  internal  charging 
calculation  for  each  scenario  and  compute  statistics  from  its  results.  Nonetheless,  for  comparison  to 
other  models,  it  is  necessary  to  present  the  worst  case  computed  as  a  running  average  flux.  Included  in 
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Figure  12:  Comparison  of  TEM-2  scenarios  and  GOES-10  observed  statistical  distribution  of  >  2  MeV 
electron  flux.  The  plot  format  follows  that  of  Figure  1.  GOES- 10  flux  is  multiplied  by  37t  as  in  Figure  11. 
The  TEM-2  distribution  shown  is  taken  from  200  Monte  Carlo  scenarios. 
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the  figure  are  GEO-specific  models,  based  on  observations  at  GEO:  IGE  [Sicard-Piet  et  al. ,  2008], 
Fennell  et  al  [2000],  and  O’Brien  et  al.  [2007]  worst  case  observed  and  flux  limits.  The  primary  curves 
of  interest  are  the  percentiles  from  40  Monte  Carlo  scenarios  in  TEM-2.  These  curves  represent  the 
best  estimate  (median),  and  various  confidence  values  from  the  Monte  Carlo  scenario  method.  Eor 
example,  there  is  nominally  a  5%  chance  that  the  average  flux  for  any  given  10  year  mission  will  exceed 
the  curve  labeled  “95%  Monte  Carlo.” 

In  both  panels  of  Figure  13,  for  reference,  TEM-2  curves  are  plotted  for  the  static  mean,  median 
and/or  95th  percentile.  However,  these  curves  do  not  necessarily  represent  realistic  environments  and 
are  really  only  useful  as  diagnostics.  Additionally,  there  is  the  comparable  result  at  GEO  from  the 
l'EM-2  reanalysis  (see  sections  3-5).  The  appropriate  static  CRRESELE  model  is  also  provided  in  each 
panel  for  comparison. 

Finally,  as  shown  in  the  top  panel  of  Figure  13,  for  a  “quick-and-dirty”  calculation  of  the  average  flux, 
one  can  avoid  the  time-dependent  scenario  calculation,  and  fly  the  mission  through  a  static  “Perturbed 
Mean"  for  each  scenario.  This  calculation  is  expedient  because  it  only  requires  the  perturbation  of  the 
flux  map  (41)  for  each  scenario  and  a  relatively  small  number  of  orbits  through  the  average  from  that 
perturbed  map,  rather  than  performing  the  full  mission-length  calculation  required  for  the  Monte 
Carlo  curves.  However,  the  “Perturbed  Mean”  neglects  the  “weather”  dynamics  (38),  and  can  give 
incorrect  results  for  short  missions.  Further,  the  “Perturbed  Mean”  can  not  be  used  for  quantities  that 
are  time- sequence  dependent  (e.g.,  a  worst  case  24-hour  flux). 

Several  points  are  evident  from  Figure  13.  First,  it  is  obvious  that  FEM-2  is  quantitatively  different 
from  AE8  and  the  other  models.  Given  that  TEM-2  used  no  GEO  data,  the  obvious  conclusion  is  that 
TEM-2  is  likely  inferior,  at  least  to  IGE  and  the  Fennell/O’Brien  models.  We  assume  that  in  future 
models  (i.e.,  AE9),  GEO  data  will  be  included  and  agreement  with  these  models  will  improve.  Second, 
the  spread  from  the  median  to  the  95th  percentile  is  only  a  factor  of  2-4,  which  is  consistent  with  the 
spread  among  the  other  models,  where  provided.  However,  the  models  do  not  always  agree  with  each 
other  within  these  error  bars.  Later  models  (i.e.,  AE9)  promise  improved  treatment  of  the 
measurement  errors,  both  including  error  estimates  for  each  measurement  and  including  more  data 
sets.  Finally,  the  worst  case  estimates,  although  in  major  disagreement  with  the  Fennell  and  O’Brien 
models  [Fennell  et  al. ,  2000;  O’Brien  et  al. ,  2007],  are  nonetheless  bounded  by  these  models,  except  at 
very  high  energies  (>  4  MeV). 


2.3.2  GPS  example 

Figure  14  shows  the  same  quantities  as  in  Figure  13,  except  for  a  nominal  GPS  orbit.  In  place  of  the 
IGE  model,  we  have  used  the  MEO  GNSS  model,  provided  by  ONER  A  [Sicard-Piet  et  a/.,  2006; 
Bourdarie  et  a/.,  2010].  For  the  average  environment  (top),  quantitative  agreement  is  surprisingly 
good.  For  the  worst  case  environment  (bottom),  the  only  comparison  is  to  CRRESELE  Max.  So,  we 
cannot  learn  much  about  quantitative  performance,  but  the  local  maxima  evident  in  the  percentiles  of 
the  Monte  Carlo  scenarios  are  very  suspicious,  and  likely  erroneous.  Since  they  are  absent  from  the 
static  95%  environment,  they  are  likely  an  artifact  of  either  large  measurement  errors  on  specific 
sensors,  or  truncation  of  the  principal  components  in  the  prior  model. 
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TEM-2,  GEO  -  Mean,  10  years 


10 


10 


10 


1 - !lll 

AE8Min 
AE8Max 
CRRESELE  Avg 
IGE  low 
IGE  mid 
IGE  high 


10 


| 

2>  104 


10 


II 


551 


s 


'  —  -  Static  Mean  Env 

■  —  -  Static  Median  Env 

■  —  -  Static  95%  Env 

■  "Reanalysis 
“ ■  50%  Monte  Carlo 
" —  75%  Monte  Carlo 
— 90%  Monte  Carlo 
“ "  95%  Monte  Carlo 
"O—  50%  Perturbed  Mean 
•O"  75%  Perturbed  Mean 
•O—  90%  Perturbed  Mean 
•<*•••  95%  Perturbed  Mean 


n 


> 

<u 

2 


(/> 


E 


o 

* 


MeV 


Figure  13:  10-year  environments  at  GEO.  lop:  Mean  flux  from  various  models;  AE8  (MIN  and  MAX 
are  the  same  at  GEO).  Bottom:  worst-case  24-hour  exponentially- weighted  flux  average. 
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TEM-2,  GPS  -  Mean,  10  years 


Figure  14:  10-year  environments  at  GPS.  Top:  Mean  flux  from  various  models;  AE8.  Bottom:  worst-case 
24-hour  exponentially- weighted  flux  average. 
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2.3.3  Dose  depth  curves 


Finally,  one  of  the  most  widely-requested  quantities  that  current  models  cannot  produce  for  an 
arbitrary  orbit  is  “error  bars  on  dose-depth  curves.”  Figure  15  show's  just  this  quantity  from  TEM-2 
for  GEO  (top)  and  GPS  (bottom)  orbits.  For  comparison,  we  have  also  plotted  dose  computed  from 
reference  mean  environment  models  used  in  Figures  13  and  14.  Again,  quantitative  agreement  is  poor, 
especially  beyond  100  mils  of  Aluminum  shielding.  Nonetheless,  TEM-2,  along  with  TEM-1,  represents 
the  first  class  of  global  trapped  radiation  model  that  can  compute  dose  “error  bars”  which  have  real 
probabilities  associated  with  them,  accounting  for  both  environmental  dynamics  and  measurement 
error.  This  is  a  major  step  toward  truly  usable  “next  generation  radiation  specifications.”  However, 
much  work  obviously  remains  to  be  done  before  a  model  like  TEM-2  can  be  used  quantitatively  for 
satellite  design. 
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TEM-2,  GEO  -  Dose,  3650  days 
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Figure  15:  10-year  total  dose  environments  at  GEO  (top)  and  GPS  (bottom). 
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3  Reanalysis  procedures 


In  addition  to  the  engineering  application  demonstrated  above,  TEM-2  can  also  be  used  as  the 
“background”  model  for  a  long-term  reanalysis.  We  will  discuss  how  the  prior  model  is  used  with 
long-term  observations  of  typically  lower  resolution  to  estimate  the  global  state  of  the  electron 
radiation  belts  in  6-hour  snapshots  over  an  entire  solar  cycle-a  statistical  reanalysis.  Whereas  data 
assimilation  usually  involves  a  maximum  likelihood  state  that  balances  observational  errors  with  the 
output  of  a  global  simulation,  in  the  following  demonstration,  the  statistical  prior  model  from  TEM-2 
will  stand  in  for  the  output  of  the  physical  simulation.  We  will  also  discuss  the  “cross-calibration”  of 
the  reanalysis  data  sets,  and  the  validation  of  the  reanalysis.  Finally,  we  will  discuss  some  potential 
scientific  uses  of  our  reanalysis. 


3.1  Stating  the  data  assimilation  problem 


The  objective  of  our  data  assimilation  is  to  obtain  the  most  likely  state  of  the  electron  belts  given  a  set 
of  observations  taken  near  a  fiducial  time  of  interest.  In  an  ideal  situation,  one  would  have  a  large  set 
of  observations  y  of  the  phenomenon  of  interest  x  (i.e.,  radiation  belt  fluxes  on  a  global  grid),  and  one 
could  solve  for  x  by  least  squares  or  more  sophisticated  inversion  methods.  Formally, 

y  =  H(x),  (42) 

where  //j(x)  is  the  “measurement  function”  for  the  ith  observation,  ?/, .  In  the  case  of  the  radiation 
belts,  the  observations  are  typically  integrals  over  flux  (or  phase-space  density)  plus  a  background  (6). 
Since  integrals  are  linear  operators,  we  can  represent  (42)  as 

y  =  b  4-  Hx.  (43) 


A  least-squares  solution  is  tempting: 

X  =  (ILtIL)  - 1  Kt  (y  -  b) .  (44) 

However,  H_  is  typically  underdetermined  (there  are  more  x’s  than  y' s),  so  that  Hj £L  is  singular.  A 
somewhat  more  sophisticated  solution  employs  the  pseudo- inverse: 

x  =  Kf(y-b).  (45) 

This  solution  at  least  exists;  however  it  solves  the  underconstrained  problem  by  minimizing  the  length 
of  the  unspecified  portion  of  x.  That  is,  it  assumes  that  where  there  is  no  measurement,  there  is  no 
flux  (or  phase-space  density).  Obviously,  this  is  not  a  physically  reliable  assumption.  Something  more 
is  needed  to  estimate  x. 

In  the  inversion  of  Energetic  Neutral  Atom  (ENA)  images,  it  is  common  to  enforce  a  smoothing 
criterion  [DeMajistre  et  a/.,  2004:  Perez  et  a/.,  2004].  Smoothing  appears  to  work  well  for  ENAs,  but 
probably  would  not  work  well  in  the  radiation  belt~our  underlying  measurements  are  fluxes  along 
spacecraft  trajectories,  rather  than  global  images.  Instead,  we  employ  a  maximum-likelihood  approach: 
we  want  to  estimate  the  most  likely  instantaneous  state  of  the  electron  radiation  belts  consistent  with 
not  only  a  set  of  observations  but  also  a  prior  model  of  the  likelihood  of  the  various  states  of  the  belts. 

We  have  elected  to  work  with  the  natural  logarithm  of  flux  (log.?)  because  flux  errors  are  typically 
treated  as  multiplicative  (e.g.,  the  traditional  “factor  of  2”)  rather  than  additive.  Also,  maximum 
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likelihood  log  flux  provides  more  realistic  solutions  than  maximum  likelihood  flux,  because  the  mode 
(i.e.,  maximum  likelihood  value)  of  most  flux  distributions  tends  toward  zero,  and  this  is  especially 
problematic  in  a  multivariate  flux  distribution.  The  mode  of  the  log  flux  distribution  is  close  to  its 
median,  thereby  representing  a  more  “typical”  value. 

3.2  The  likelihood  to  be  maximized 


Accepting  that  we  will  compute  the  maximum  likelihood  of  log x  constrained  to  match  a  set  of 
observations  y ,  in  the  usual  notation  of  Bayesian  conditional  probabilities  we  want  to  find  the  value  x 
that  maximizes: 

£(logf)  =  p(log  x\y)  =  (46) 

p{y) 

Unfortunately,  we  do  not  usually  know  p(logx, y)  or  p(y).  However,  there  are  two  quantities  that  we 
can  obtain,  from  which  we  can  compute  enough  of  p(log x\y)  to  find  a  maximum  likelihood  solution. 
Namely,  we  know  the  part  of  £(logx)  that  depends  on  x: 


p(logx|jf)  = 


oc 


p(\ap,.r.y) 

p(y) 

p(v I  logx)p(logf). 


(47) 


In  the  last  expression,  we  havre  dropped  th ep(y)  factor  because  it  is  invariant  with  respect  to  x.  1  he 
two  quantities  that  remain  are  p(y\  logx),  which  is  the  probability  of  a  given  set  of  measurements  (y) 
assuming  we  know  the  state  of  the  electron  belts  (x),  and  p(logx),  which  is  onr  so-called  “prior”  model. 
We  note  that  p(y\  logx)  =  p(y|x)  because  the  condition  on  logx  is  equivalent  to  the  condition  on  x. 


As  is  usually  the  ca.se,  instead  of  actually  maximizing  £(logx),  we  will  minimize  its  negative  logarithm 
(without  the  terms  that  are  constant  with  respect  to  x): 

-  l'ig^(logx)  ~  £(logx)  =  -  log p{y\  logx)  -  logp(logx).  (48) 


The  error  covariance  matrix  for  the  maximum  likelihood  logx  is  approximated  by 


/ 

il  =  COY  log  X  = 

V 


d2£ 

d  log  xxd  log  Xj 


Idle  standard  error  ^Xi  on  Xi  is  given  by: 


) 


(49) 


(50) 


3.3  The  statistical  prior  model  as  a  likelihood  function 


We  begin  to  construct  the  terms  in  the  maximum  likelihood  problem  (48)  with  the  “prior”  model 
p(logx).  We  will  start  with  p(q)  and  work  our  way  to  p(logx).  We  know  from  section  2.2.8,  that  <f  lias 
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a  multivariate  Gaussian  distribution: 


p(q)  = 


exp[-(f'  g/2] 

( 2ir)N« /2 


(51) 


In  the  following  we  will  need  the  generic  relationship  that  describes  a  multivariate  change  of  variables 
in  probability  calculus: 


p(v)  =  p(w) 


dw 

dv 


(52) 


where  \dw/dv\  is  the  determinant  of  the  Jacobian  of  the  change  of  variables.  Further,  we  will  need  the 
related  expression  for  monotonically  increasing  univariate  changes  of  variables 


dw  p(v) 
dv  p(u’) 


(53) 


Both  of  these  expressions  can  be  derived  from  the  conservation  of  probability  density,  which  is 
analogous  to  the  conservation  law  that  governs  changes  of  variables  in  integrals  or  phase  spaces. 


Applying  the  conversion  from  q  to  z ,  we  have: 

../-V  _  exp[-zI'|T1z/2] 
P\Z)  r - -  J 

v/|S|(2ff)^/2 

noting  that  |E|  =  \Q\2  in  the  idealized  case  of  a  non-singular  E. 


(54) 


Converting  to  x,  we  have: 


exp  [-5^5  1?/2]rr/i(xi) 


(55) 


where  fi  is  the  probability  density  function  for  x*,  and  <\>  is  the  probability  density  function  for  a  unit 
normal  (standard  Gaussian): 


/<(*«) 


0(z) 


Oxi  ’ 

exp  [— z2/2]  d<i> 

v/27t  dz ' 


(56) 

(57) 


Finally,  converting  to  log x,  we  have: 


„  ^  exp[-zr^lz/2 ]1-rxj/i(xi) 

p(log  j)  =  — - —  ~ -  II  ,,  \ 

^f(27r)‘v-/2  i  z *) 

When  we  take  the  negative  logarithm  and  replace  ^E-1  z  with  qT  q,  we  have 

-  log p(log  x)  =  -T:  [log  X,  +  log  /, (x,)] 

i 

+ i  {qrq  —  z1'  z)  -I-  constants. 


(58) 


(59) 


We  note  that  the  transform  from  z  to  q  is  constant  and  linear,  and  therefore  it  does  not  modify  the 
log-likelihood  function,  except  by  the  addition  of  a  constant,  which  drops  out  of  the  optimization. 
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(60) 


For  fi  following  a  Weibull  distribution,  we  differentiate  (11) 

c<  \ 

/(x;x0,7)  =  — —exp 

*0 

which  leads  to 


-logp(logx) 


-  [7i  log 27  -  (x»/xo,iP'l 

i 

+  -  (q1  q  —  zr z)  +  constant s. 


For  fi  following  a  log-normal  distribution  (14),  we  have  a  remarkably  simple  result: 

logJ'j  - 


Ozi 

OjCi 

fi(xi) 


-  logp(Iogx) 


(XjfT,) 

1  ^  /  log  Xj  -  pj\  _  <A(gj) 

xi°i  \  n\  ) 

q  4  constants. 


(61) 


(62) 

(63) 

(64) 


(65) 


Depending  on  our  choice  of  marginal  distribution,  either  ((>1)  or  (65)  supplies  —  logp(logx)  in  the 
maximum  likelihood  problem  (48).  In  either  case,  the  vector  of  medians  occurs  at  2  =  0,  which  is  also 
q  —  0,  providing  a  reasonable  starting  point  for  the  numerical  optimization  algorithm.  The 
optimization  will  occur  in  q \  but  will  maximize  the  likelihood  of  log x. 

In  the  reanalysis,  we  use  a  lognormal  marginal  distribution  rather  than  the  Weibull  used  in  the  rest  of 
l  lv\I-2.  1  lie  Weibull  tends  to  pull  the  reanalysis  solution  away  from  the  data  because  the  maximum 
likelihood  log  flux  at  some  gridpoints  is  much  higher  than  the  median  log  flux.  By  using  a  lognormal, 
for  which  the  maximum  likelihood  log  flux  is  tin'  median  flux  we  effectively  penalize  deviations  from 
the  median,  rather  than  from  the  maximum  likelihood  log  flux. 


3.4  The  grid  interpolation  functions 


In  order  to  relate  the  radiation  belt  state  (q  or  .r)  to  observations  (i.e.,  to  compute  p(y\  logx),  we  will 
need  to  interpolate  between  grid  points.  1  he  flux  x  at  any  given  E,  acq ,  and  Lrn  is  given  by  the  sum  of 
a  set  of  tri-linear  finite  elements: 


x(Ij  ,  Orq,  Lm)  —  ^  XjVi(E)  Ctegi  f-'m)* 


Each  element  is  the  product  of  three  linear  contributions: 

Vi(E,acq,  Lrn)  =  v\E\E)v\a\aeq)vlL)(Lm), 


(66) 


(67) 


and  the  linear  contributions  are  defined  as: 


„(*) 


(X) 


x-xr 

vhiKh  \* 
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mid 
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rhiRh  _ 
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<  X  <  X 


high 


ot  herwise 


(68) 
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for  X  being  either  E ,  ae<?,  or  Lrn. 

The  mid-point  of  one  element  is  the  high  or  low  point  of  its  neighbor,  as  would  be  expected.  At  the 
atmospheric  boundary,  the  low  point  is  set  to  the  Lm-dependent  analytical  loss  cone  angle  provided  by 
(1),  and  no  grid  points  exist  inside  the  loss  cone.  The  flux  is  zero  beyond  the  mid-points  of  the  grid 
points  at  the  other  boundaries  of  the  grid. 


3.5  Measurement  functions 


We  consider  a  measurement  to  be  an  integration  over  time  and  velocity  space.  That  integration  will  be 
discretized  on  our  grid,  and  represented  as  a  weighted  sum  over  the  model  fluxes  plus  an  expected 
background  (6).  Writing  equation  (43)  in  subscript  notation,  we  have: 

yk^bk+Y^,  HkiXi-  (69) 


As  we  w  ill  see,  the  background  can  be  estimated  from  other  sensors  on  the  vehicle.  For  the  calculation 
of  //,  we  ignore  temporal  dependence  over  the  duration  of  the  measurement  (7*),  and  we  represent 
velocity  space  in  energy  E  and  local  pitch-angle  a  coordinates.  Thus,  the  measurement  function  for 
observation  y (in  counts)  taken  at  (Lm,  B/Beq)  is  given  by: 
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where  aeq(a)  relates  equatorial  and  local  pitch-angle  via  conservation  of  the  first  adiabatic  invariant: 

(71) 
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The  bounds  of  pitch-angle  integration  a\  and  Q2  are  local  pitch  angles  that  correspond  to  either  the 
low  and  high  equatorial  p>itch  angle  limits  of  the  grid  point,  or  a2  =  tt/2,  when  the  high  equatorial 
pitch  angle  limit  of  the  grid  point  includes  locally  mirroring  particles.  The  function  Ak(E,  a)  provides 
the  instrument’s  effective  area  (in  cm2)  for  particles  with  energy  E  and  local  pitch-angle  a.  The  units 
of  II ki  are  (cm2srskeV).  We  have  chosen  to  represent  yk  as  counts  rather  than  counts- per-second  or 
flux  because  the  number  of  counts  is  a  more  accurate  depiction  of  the  usual  measurement  process,  and 
it  allows  us  to  approximate  the  measurement  error  through  a  Poisson  counting  process  when  there  are 
fewr  counts.  In  our  maximum  likelihood  treatment,  we  will  use  two  error  distributions,  depending  on 
w'hether  Poisson  counting  error  dominates  instrument  resp>onse  uncertainty. 


3.6  Pseudo-Poisson  error  distribution 


If  wre  assume  that  the  kth  measurement  is  a  Poisson  process  w  ith  expected  counts  A*,  we  have: 
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(72) 

(73) 


34 


We  next  introduce  /?fc,  a  variance  inflation  factor  (see  section  4.3),  which  will  allow  us  to  maintain  the 
shape  of  the  Poisson  distribution  while  modifying  its  variance.  We  will  do  this  when  we  must  optimize 
(48)  for  a  variety  of  data  points  that  are  not  all  taken  simultaneously.  We  can  preserve  the  maximum 
likelihood  value  of  p(yk\^k)  while  inflating  the  variance  by  1  /j3k  using  the  following  transform: 

p(yk\^k'<Pk)  =  'PiykPkl^kPk)-  (74) 


1  he  negative  log  and  gradient  are  given  by: 

1=  —  logp(y|A;/?) 

=  j3(X  —  y  log  A)  4-  constants, 

(75) 
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3.7  Gaussian  error  distribution 


For  most  of  our  data,  for  sufficiently  large  counts,  the  predominant  error  source  is  our  lack  of  knowledge 
of  the  true  measurement  function,  II  (x).  It  is  common  to  represent  that  uncertainty  as  relative  error 
(<fy),  implying  an  approximately  Gaussian  distribution  of  the  error  in  log  counts  or  count  rate. 


I  he  Gaussian  error  distribution  for  the  natural  log  of  yk  is: 

yk 
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Again,  we  introduce  a  variance  inflation  factor  3^,  which  appears  simply  as 

£(logyfc|  log  Xk;6yk/y/%) 


yk 


(78) 


The  negative  log  and  gradient  are  given  by: 
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Table  3:  Data  channels  used  to  drive  the  reanalysis 
Channel  Sensor  Type  Energy  6y 


HEO-1,  1994-2008,  L  >  4 

El 

Telescope 

>0.13  MeV 

1.4 

E2 

Telescope 

>0.23  MeV 

1.2 

E3 

Dosimeter 

>1.5  MeV 

1.2 

E4 

Dosimeter 

>4.0  MeV 

1.1 

E5 

Dosimeter 

>6.5  MeV 

1.4 

E6 

Dosimeter 

>8.5  MeV 

1.1 

HEO-3,  1997-2008,  L  >  2 

E3 

Telescope 

>0.45  MeV 

1.5 

E5 

Dosimeter 

>1.5  MeV 

1.3 

E6 

Dosimeter 

>3.0  MeV 

1.3 

ICO,  2001-2008,  L  >  2.5 

El 

Dosimeter 

>0.95  MeV 

0.93 

E2 

Dosimeter 

>2.0  MeV 

1.0 

E3 

Dosimeter 

>3.5  MeV 

0.97 

E4 

Dosimeter 

>5.5  MeV 

0.79 

E5 

Dosimeter 

>6.8  MeV 

0.81 

SAMPEX,  1992-2004,  all  L 

PET/ELO  Telescope 

2-6  MeV 

* 

*  see  section  4.2 


4  Reanalysis  data 


For  the  reanalysis,  we  rely  primarily  on  a  combination  of  simple  sensors  developed  by  The  Aerospace 
Corporation  and  known  as  IIEO-1,  HEO-3,  and  ICO.  Each  of  these  sensors  can  be  thought  of  as  an 
approximately  omnidirectional  integral  electron  channel.  The  HEO  vehicles  are  in  12-hour,  critically 
inclined,  Molniya  orbits,  giving  them  access  to  Lm  ~  2  at  the  magnetic  equator  and  higher  L  shells  at 
higher  magnetic  latitudes.  The  IIEO-1  vehicle  only  telemeters  data  for  Lm  >~  4.  The  ICO  vehicle 
(2001-026)  is  in  a  circular  orbit,  at  ~ 10,390  km,  with  an  inclination  of  45°,  crossing  the  magnetic 
equator  at  Lm  ~  2.5. 

Finally,  to  provide  some  indication  of  the  variation  of  inner  zone  electrons,  we  utilize  the  SAMPEX 
PET/ELO  “L  bin”  history  provided  by  S.  Kanekal  (private  communication),  from  1992-2004. 
SAMPEX  is  in  a  circular  ~  600  km  orbit  at  82°  inclination.  In  section  (4.2),  we  will  discuss  how  this 
satellite  that  is  primarily  in  the  loss  cone  is  used  to  constrain  the  trapped  population. 


4.1  HEO  and  ICO  data 

Both  HEO-1  and  HEO-3  carry  simple,  two-element  solid-state  telescopes  (El  and  E2  in  Table  3). 
However,  on  HEO-3,  these  channels  are  unusable  for  years  at  a  time  due  to  what  appears  to  be  a 
thermal  background.  Further,  HEO-3  carries  an  additional  telescope,  E3,  with  minimal  collimation. 
The  attitude  information  for  these  vehicles  is  not  available,  so  the  telescopes  are  treated  as 
omnidirectional  sensors.  Both  HEO  vehicles  and  ICO  carry  dosimeters,  consisting  of  a  silicone  target 
behind  a  dome  shield.  Each  dosimeter  has  a  unique  dome,  to  distinguish  incident  particle  energy,  and 
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Figure  16:  Calculated  energy  response  curves  for  HEO  and  ICO  channels. 
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analysis  of  the  energy  deposited  in  the  silicon  target  provides  rudimentary  differentiation  between 
protons  and  electrons.  For  each  “electron”  channel,  there  is  a  corresponding  “proton”  channel  that  is 
used  to  estimate  the  channel’s  background.  Before  being  used  in  the  reanalysis,  the  HEO  and  ICO 
data  are  binned  into  0.25  Lm  bins  using  the  OPQ  magnetic  field  model. 

Although  we  have  not  formally  intercalibrated  the  HEO  and  ICO  data,  we  did  perform  some 
preliminary  corrections  based  on  the  observed  response  to  galactic  cosmic  rays  when  the  vehicles  are 
outside  the  radiation  belts.  This  led  to  some  constraint  on  the  relative  geometric  factors  of  the  sensors. 
Further,  we  used  early  iterations  of  the  reanalysis  to  correct  for  some  apparent  systematic  errors  in  the 
channel  response.  For  example,  if,  in  sample,  the  reanalysis  solution  was  systematically  lowrer  than  the 
observed  counts  in  a  particular  channel,  wre  assumed  that  channel’s  response  function  needed  to  be 
rescaled  upward. 

As  for  the  background,  a  simple  scatter  plot  of  each  dosimeter  electron  channel  to  its  corresponding 
proton  channel  (not  showrn)  reveals  a  linear  cutoff  for  lowest  electron  count  rate  as  a  function  of  the 
proton  count  rate.  By  assuming  (safely)  that  this  cutoff  represented  the  case  wrhere  the  proton 
background  dominated  the  electron  signal,  we  were  able  to  determine  a  simple  scale  factor  (a  multiplier 
of  ~  2  —  3)  converted  the  count  rate  in  each  proton  channel  into  a  proxy  background  for  its 
corresponding  electron  channel. 


4.2  SAMPEX  data 


Data  from  SAMPEX  is  somewdiat  problematic  because  it  is  nearly  always  in  the  loss  cone.  However, 
large-scale  global  correlations  [Kanekal  et  ai ,  2001]  in  the  electron  belts  suggest  that  such  a  sensor  can 
nonetheless  provide  valuable  information  about  the  gross  behavior  of  the  otherwise- unobserved  inner 
belt  electrons.  To  utilize  the  SAMPEX  data  we  relate  it  to  the  HE03/E5  (>  1.5  MeV)  channel  for  the 
high  altitude  legs  of  the  HEO-3  orbit.  For  each  0.1  Lm  from  3  to  8,  we  obtain  powrer-law  regression 
parameters  between  SAMPEX  PET/ELO  [Cook  et  al .,  1993]  and  HE03/E5  based  on  daily  averages. 

In  addition  to  the  power-law  parameters,  we  compute  the  typical  B/Beq  value  of  the  HE03/E5  data, 
and  we  compute  an  RMS  error  in  the  log  fluxes,  which  serves  as  Sy.  Figure  17  shows  the  Lm 
dependence  of  these  parameters.  For  Lm  <  3,  where  the  powrer-law  regression  breaks  down,  wre  use  the 
regression  values  for  Lm  =  3,  but  multiply  Sy  by  3,  to  produce  only  qualitative  influence  over  the  inner 
zone.  The  background  is  assumed  to  be  0,  and  a  Gaussian  relative  error  distribution  is  used. 


4.3  Accounting  for  latent  data 


In  order  to  make  the  reanalysis  robust  to  data  gaps,  we  include  delta  outside  the  nominal  time  window 
for  a  particular  fiducial  time.  That  is,  the  fiducial  cadence  is  6  hours,  but  we  include  data  from  a  larger 
window.  The  minimum  time  window  includes  2  whole  passes  through  the  radiation  belt.  However,  in 
case  of  a  data  gap,  the  window  is  expanded  in  either  or  both  directions  in  time  to  allowr  at  least  one 
full  pass  before  and  one  after  the  fiducial  time. 

To  degrade  the  influence  of  data  that  is  farther  in  the  past  or  future,  wre  need  to  knowr  how  the 
autocorrelation  scales  with  time  offset.  In  Figure  18,  wre  have  calculated  the  rank  order  correlation 
coefficient  (Spearman’s  p)  between  electron  fluxes  in  each  HEO-3  channel  and  their  counterparts  1  day 
later  (same  orbit  leg  and  odd /even  parity).  We  see  that  there  is  very  little  dependence  on  energy  or 
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Figure  17:  (Top)  SAMPEX  PET/ELO  intensity  can  be  used  as  a  proxy  for  HEO-3  E5.  (Bottom)  The 
parameters  of  the  proxy  relationship  vary  with  L  shell. 
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Table  4:  Coefficients  for  1-day  autocorrelation  function  vs  L  for  HEO-3  fluxes 


^24,1 
T  24,2 
7*24,3 
7*24,4 


0.504938 

0.477943 

1.130208 

4.254812 


orbit  status,  at  least  at  the  higher  energies.  The  correlation  coefficient  can  be  fit  using  the  function: 


T2i  (L)  =  r2  4,1  +  r24j2  exp 


(81) 


where  the  coefficients  are  given  in  Table  4. 

Figure  19  shows  that  the  correlation  at  2  days  is  well  modeled  by  the  square  of  the  correlation  at  1  day. 
This  is  consistent  with  a  first-order  autoregressive  process  with  correlation  at  arbitrary  time  difference 
At  (in  days)  given  by: 


7*heo (L,At)  ~  r24(T)*A^ 


(82) 


With  a  definition  of  how  the  autocorrelation  varies  with  time  offset,  we  can  define  the  variance 
inflation  factor  j3 k  for  measurement  y 


Pk  —  ruEo(^ki  Atk)-  (83) 

This  (3k  is  used  in  the  measurement  penalty  functions  given  in  sections  3.6  and  3.7.  Comparing  to 
equations  (75)  and  (79),  we  see  that  /3  is  a  linear  multiplier  on  the  likelihood  function.  Therefore,  as 
the  correlation  goes  to  zero  with  increasing  time  offset,  the  effect  of  the  measurement  in  the  overall 
likelihood  function  goes  to  zero  also. 
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R(24)  =  0.50  +  0.48  exp(-[(L-1 ,13)/4.25]4) 


HEO-3  Electrons  Autocorrelation 


Orbit  Status: 

0:  Even,  out,  high 
1:  Odd,  out,  high 
2:  Even,  in,  high 
■  3:  Odd,  in,  high 
4:  Even,  out,  low 
5:  Odd,  out,  low 
6:  Even,  in,  low 
7:  Odd,  in,  low 


'ignre  18:  Autocorrelation  for  HEO-3  electron  fluxes  at  1  day  lag  for  all  orbit  status  and  energy  channels. 
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Figure  19:  Autocorrelation  for  HEO-3  electron  fluxes  at  2  days  lag  versus  the  square  of  the  correlation 
at  1  day  lag  for  all  orbit  status  and  energy  channels.  The  good  agreement  suggests  that  the  HEO-3 
autocorrelation  is  a  first-order  autoregressive  process. 
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Channel 

<  logio 

Table  5:  ICO  out-of-sample  validation  statistics 
(y/A)  >  <  log10  |y/A|  >  RMSE(log10  y,  log10  A) 

PE  % 

std  logio  y 

El 

0.410 

0.572 

0.771 

74.1 

1.515 

E2 

0.490 

0.628 

0.816 

72.5 

1.556 

E3 

0.297 

0.633 

0.789 

69.0 

1.417 

E4 

0.048 

0.479 

0.612 

66.4 

1.055 

E5 

0.071 

0.339 

0.436 

73.8 

0.853 

5  Reanalysis  results 


We  performed  the  reanalysis  process  for  every  6  hour  time  step  from  1992-2008.  Figure  20  shows  the 
entire  interval,  and  two  shorter  sub- intervals.  One  can  see  that  there  are  dynamics  in  all  three  regions 
of  the  electron  belts:  inner  belt  evolution  (driven  by  SAMPEX  data),  slot  filling  and  emptying,  and 
outer  zone  dynamics.  The  progressive  filling  and  emptying  of  the  slot,  in  particular,  is  striking.  One 
can  also  observe  several  occasions  when  the  inner  edge  of  the  outer  zone  retreats  to  higher  L  during 
magnetically  quiet  periods.  These  phenomena  are  fairly  well  known,  so  their  presence  is  reassuring. 


5.1  Validation  and  known  issues 


To  validate  the  reanalysis  results,  we  have  performed  the  reanalysis  procedure  on  a  shorter  interval 
(2002  and  2003)  without  ICO  data.  We  can  then  reconstruct  what  ICO  would  have  seen  in  that 
reanalysis  and  compare  with  observations.  Figure  21  shows  several  comparisons  of  ICO  observations 
( y )  to  the  reanalysis  (A)  in  counts.  In  the  top  left  panel,  we  see  in-sample  validation:  when  ICO  data  is 
used  in  the  reanalysis,  how  well  does  the  optimized  flux  map  reproduce  the  ICO  data?  We  see  that 
even  in  sample,  the  presence  of  other  data  sets  pulls  the  solution  away  from  the  ICO  data  to  some 
extent,  especially  at  higher  Lm  where  the  magnetic  coordinates  are  least  reliable,  the  other  panels  of 
Figure  21  show  out-of-sample  comparisons  for  each  ICO  data  channel.  The  lower-energy  channels  (El, 
E2)  show  noticeable  overestimates  in  the  reanalysis  for  Lm  ~  5  —  6.  One  explanation  for  this 
disagreement  could  be  the  reduced  performance  at  high  Lm  of  the  magnetic  coordinates  we  have 
chosen,  and  another  explanation  could  be  incomplete  intercalibration  of  the  IIEO  and  ICO  data. 

Quantitatively,  we  can  compare  the  observed  ( y )  and  estimated  (A)  ICO  counts  in  a  number  of  ways, 
table  5  provides  out-of-sample  validation  statistics  for  each  ICO  channel.  The  <>  operator  defines  a 
sample  average,  and  PE  stands  for  prediction  efficiency,  which  is  the  ratio  of  the  mean  squared  error 
between  log 10y  and  log10  A  to  the  variance  of  log10y  subtracted  from  one,  expressed  as  a  percentage. 

1' he  final  column  provides,  for  reference,  std  log10  y,  the  standard  deviation  of  log10  y.  1  he 
performance,  in  terms  of  prediction  efficiency,  is  comparable  across  all  channels.  1  he  mean  absolute 
error  (column  3)  and  the  root  mean  squared  error  (RAISE)  (column  4)  correspond  to  multiplicative 
errors  on  the  order  of  a  factor  of  2.  The  fact  that  the  mean  error  in  column  2  is  not  zero  suggests  that 
A  is  systematically  higher  than  y,  by  30%-50%  for  E1-E3,  but  by  less  than  10%  for  E4  and  E5.  These 
systematic  errors  suggest  incomplete  inter  calibration. 

A  more  detailed  survey  of  other  energies  and  pitch  angles,  besides  what  is  shown  in  Figure  20,  reveals 
some  expected  shortcomings  of  the  reanalysis.  1  he  reanalysis  exhibits  poor  performance  when  no 
bounding  data  are  available;  e.g.,  in  the  inner  zone  away  from  1.5  MeV  or  in  the  outer  zone  before 
HEO-3  data  begins.  Further,  since  the  constraint  data  are  nominally  omnidirectional  integral  channels, 
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Figure  20:  Flux  versus  Lm  profiles  versus  time  in  the  TEM-2  Reanalysis:  (top)  the  entire  reanalysis, 
(middle)  one  year,  centered  on  the  July  2004  storm,  and  (bottom)  one  month  centered  on  July  24,  2004. 
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Figure  21:  Validation  of  TEM-2  reanalysis  versus  ICO  (in  and  out  of  sample). 
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the  reanalysis  sometimes  produces  unphysical  or  at  least  suspicious  differential  energy  spectra  and 
angular  distributions.  This  problem  likely  stems  primarily  from  a  lack  of  intercalibration,  but  it  also  is 
compounded  by  our  choice  of  magnetic  coordinates  (i.e.,  bounce  rather  than  drift  invariants  and  a 
quiet  rather  than  active  magnetic  field  model). 


5.2  Engineering  application 


With  appropriate  caution,  we  can  move  on  to  apply  the  reanalysis  results  to  engineering  applications. 
Figures  11,  13,  14,  and  15  all  include  a  trace  for  the  TEM-2  reanalysis.  While  the  reanalysis  shows 
generally  poor  quantitative  performance  (possibly  due  only  to  inter  calibration),  it  often  shows 
qualitative  agreement.  In  particular,  the  reanalysis  in  Figure  11  captures  some  of  the  sharp  increase  in 
flux  at  GOES  on  March  18-20,  2006.  The  dynamics  in  the  reanalysis  represent  real  events  with 
more-or-less  complete  temporal  correlations  (rather  than  the  first  order  temporal  behavior  represented 
statistically  in  the  Monte  Carlo  scenarios).  This  potentially  very  realistic  environmental  dynamics 
could  be  used  to  characterize  idiosyncratic  effects  of  individual  storms.  The  Monte  Carlo  scenarios  and 
the  static  models  have  no  hope  of  representing  that  type  of  idiosyncratic  behavior.  It  remains  to  be 
seen  whether  a  more  quantitatively  accurate  reanalysis  model  will  be  able  to  reproduce  the  worst  cases, 
since  the  nonlinear  optimization  (48)  tends  to  suppress  extreme  values. 


5.3  Science  application 


The  reanalysis  simplifies  the  use  of  standard  tools  developed  for  particle  measurements,  and,  as  we  will 
see,  it  enables  use  of  some  uniquely  global  tools.  The  reanalysis  allows  one  to  pick  a  particular  energy, 
equatorial  pitch  angle,  and  Lm  value  and  examine  only  the  temporal  variation  of  flux  or  phase-space 
density;  or,  one  can  select  a  slice  along  one  or  more  of  the  dimensions  to  perform  multivariate  time 
series  analysis,  without  having  to  explicitly  account  for  the  orbital  motion  of  the  observatories  or  data 
gaps  in  one  or  another  sensor  data  set.  However,  as  we  will  see,  reanalysis  also  allows  one  to  do  things 
that  were  never  before  possible. 


5.3.1  Empirical  determination  of  the  time  evolution  operator 

Space  physicists  usually  assume  that  the  phase  space  density  ft  (a  continuous  function  of  one  or  more 
adiabatic  invariants)  evolves  according  to  a  linear  operator  (usually  diffusion): 

ft=Dtft+St ,  (84) 

where  ft  is  a  time  derivative,  Dt  is  a  time-varying  linear  operator,  and  St  is  a  time-varying  combined 
additive  source  and  loss  term.  If  one  assumes  that  Dt  is  constant  at  Dj  among  a  set  of  times  T,  e.g., 
all  times  when  the  magnetic  index  Kp  is  in  some  narrow  range,  then  one  can  actually  estimate  Dj 
from  a  multivariate  sample  of  ft  for  t  6  T. 

First,  one  discretizes  ft  to  ft ,  with  time  step  At,  giving  a  matrix-vector  time  evolution  equation: 

ft+ i  —  ft  T  At  (j^ft  T St'j  , 

=  (z+  At^)  ft  +  AtSt.  (85) 
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1  lien,  for  t  £  T,  one  removes  the  average  density  and  computes  sample  correlations: 


9t  - 

(86) 

Sr  =  (ml  )r , 

(87) 

RT  —  (<jtg!+i)r  • 

(88) 

We  have  made  the  very  safe  assumption  that  ^  ft+i^  ^  • 

One  then  subtracts  \  from  both  sides  of  (85),  multiplies  the  result  by  g(  ,  and  takes  the  sample 
average  over  t  £  T ,  giving: 

( 9t+\9t  )r  =  (l  +  ^tQrr)  ( 9t9t  )r  +  ( (&tSt  -  9t  •  (89) 

Next,  we  assume  that  St  is  uncorrelated  with  /t,  i.e.,  that  additive  sources  and  losses  are  uncorrelatcd 
with  phase  space  density.  This  is  a  reasonable,  but  not  guaranteed,  assumption,  as  most  descriptions  of 
losses  in  the  literature  are  multiplicative  (e.g.,  ft/r  decay  terms),  and  most  additive  source  terms  are 
act  ually  parameterizations  of  linear  processes  operating  on  parts  of  /  outside  the  simulation  domain 
(e.g.,  a  source  due  to  stochastic  diffusion  from  lower  energies).  In  principle,  the  former  can  be 
addressed  by  incorporating  1/r  into  1),  and  the  latter  can  be  addressed  by  expanding  the  simulation 
domain,  or  by  assuming  that  the  additive  source  is  constant  over  the  times  identified  by  T  (e.g.,  the 
source  is  also  parameterized  with  Kp ,  and  would  be  represented  by  a  constant  Sq-).  1  herefore  the  last 
term  on  the  right  is  assumed  to  be  zero,  leaving: 


Hr  = 

(l  +  A/|2r)  +  0, 

(90) 

Hr  = 

(£K'-i)i' 

(91) 

In  this  formulation,  1)  need  not  be  a  diffusion  operator.  However,  if  it  is,  then  Dq.  strongly  constrains 
the  diffusion  coefficient,  which  can  be  estimated  directly  from  its  eigenvectors  [see,  e.g.,  Schulz  and 
Lanzerotti ,  1974,  section  V.3,  Spatial  Quadrature].  In  any  case,  one  can  imagine  obtaining  and  Sq- 
for  each  value  of  Kp  or  some  other  index,  and  thereby  obtaining  a  purely  empirical  model  of  the  time 
evolution  of  the  radiation  belts. 


5.3.2  Posterior  principal  component  analysis 


Another  analysis  technique  that  is  possible  only  with  the  global  time- varying  specification  provided  by 
a  reanalysis  is  principal  component  analysis.  Whereas  in  2.2.8  we  used  principal  components  (PCs)  to 
resolve  the  problem  of  the  singular  spatial  covariance  matrix,  we  will  here  define  a  set  of  posterior 
principal  components,  which  are  the  eigenvectors  of  the  spatial  covariance  matrix  of  the  log  fluxes.  As 
it  happens,  the  covariance  of  log  flux  is  invariant  to  the  choice  of  flux  or  phase  space  density.  1  he 
eigenvectors  are  not  invariant  to  the  choice  of  a  grid.  I  he  grid  we  use  for  TEM-2  (section  2) 
approximately  reflects  our  intuitive  sense  of  where  the  information  lies  in  the  radiation  belts:  it  is 
uniform  in  equatorial  pitch  angle,  logarithmic  in  energy,  and  uneven  in  Lm.  The  Lm  dimension  has 
higher  density  at  low  Lm,  and  its  midpoint  in  the  Lm  dimension  is  3.75,  roughly  at  the  inner  edge  of 
the  outer  zone;  therefore,  the  TEM-2  grid  gives  approximately  as  many  grid  points  to  the  inner  zone 
and  slot  as  it  does  to  the  outer  zone. 
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PC-1,  <xe  =90°,  36.4%  of  Variance 
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Figure  22:  The  first  posterior  principal  component  of  the  TEM-2  reanalysis,  for  equatorially  mirroring 
electrons. 


The  posterior  principal  components  we  obtain  are  individually  normalized  to  have  a  Pythagorean 
length  corresponding  to  the  number  of  variables  in  the  “reduced”  grid,  and  all  have  an  arbitrary  sign 
which  we  set  to  positively  correlate  with  solar  wind  velocity  (see  below).  The  associated  eigenvalue 
provides  the  amount  of  relative  variance  (in  log10  flux)  described  by  each  component.  We  then 
normalize  the  eigenvalues  to  sum  to  unity,  and  express  them  each  as  a  percentage. 

Figure  22  shows  the  first  posterior  principal  component  (PC-1)  for  equatorially  mirroring  electrons.  As 
is  often  the  case,  the  first  PC  has  nearly  all  entries  of  the  same  sign  and  reflects  the  global  coherence 
reported  by  Kanekal  et  al.  [2001].  In  this  case,  global  coherence  explains  36.4%  of  all  the  observed 
variance.  Global  coherence  is  likely  associated  with  the  Dst  effect  [e.g.,  Kim  and  Chan ,  1997],  an 
adiabatic  effect  not  represented  in  our  quiet  field  model  (OPQ),  and  likely  also  with  plasma- wave 
particle  interactions  that  can  span  a  large  Lm  range  simultaneously.  Specifically,  losses  associated  with 
storm- time  chorus  and  electromagnetic  ion-cyclotron  waves  can  act  on  all  Lm  values  beyond  the 
plasmapause  [Millan  and  Thome ,  2007;  Shprits  et  al .,  2008b].  During  a  magnetic  storm,  the 
plasmasphere  is  often  eroded  to  very  low  Lm  values  [e.g.,  Moldwin  et  al. ,  2003].  Inside  the 
plasmasphere,  hiss  waves  can  cause  losses  everywhere,  and  plasmaspheric  hiss  appears  to  be  driven  by 
chorus  waves  reflecting  into  the  plasmasphere  from  outside  [Bortnik  et  al .,  2008]. 

In  addition,  PC-1  shows  that  the  region  from  200  keV  to  2  MeV  near  Lm  ~  3  exhibits  a  large  degree  of 
variability.  This  reflects  the  motion  of  the  inner  edge  of  the  outer  zone,  which  is  evident  in  Figure  20. 

The  second  principal  component,  PC-2,  is  shown  in  Figure  23.  The  top  panel  is  in  the  same  format  as 
Figure  22,  and  shows  a  node  at  Lm  ~  4,  with  antinodes  at  Lm  ~  3  and  broadly  over  the  outer  zone. 
Because  the  node  is  a  horizontal  line  at  approximately  constant  Lm,  PC-2  describes  decoupling  of  the 
lower  L  shells  from  the  higher  ones.  Specifically,  this  implies  that  radial  diffusion  is  weaker  than  energy 
or  pitch-angle  diffusion,  in  a  global  context.  The  bottom  panel  illustrates  that  this  feature  is  present 
across  all  pitch  angles:  the  green  isocontour  is  the  node,  and  it  is  essentially  a  planar  surface  cutting 
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across  energy  and  pitch  angle  dimensions  at  a  nearly  constant  Lrn.  The  second  branch  of  the  node, 
folding  back  to  low  Lrn  and  progressively  higher  pitch  angles  represents  the  atmospheric  loss  cone, 
where  the  eigenvector  is  forced  to  zero  as  a  boundary.  PC-2  explains  18.2%  of  the  variance. 

Figure  24  shows  PCs  3  and  4,  which  have  somewhat  hard-to-explain  nodal  structure,  lb  some  extent  it 
is  always  tricky  to  interpret  the  higher-order  PCs  because  they  must  be  mathematically  orthogonal  to 
the  lower  order  PCs  even  though  the  principal  physical  processes  are  correlated.  In  both  PC  3  and  4 
the  node  has  more  dependence  on  energy  than  in  PC  2,  suggesting  some  decoupling  of  higher  and 
lower  energies.  PCs  3  and  4  appear  to  represent  a  combination  of  radial  and  energy  decoupling,  but 
the  interpretation  is  somewhat  ambiguous,  as  energy  and  Lm  are  not  independent  adiabatic  invariants. 
In  any  case,  there  is  no  strongly  pitch-angle  dependent  nodal  structure  in  any  of  PC-1  through  PC7, 
suggesting  that,  as  expected,  pitch  angle  diffusion  is  faster  than  both  energy  and  radial  diffusion. 

Based  on  the  PCs,  the  ordering  is:  radial  diffusion  is  slowest,  pitch  angle  diffusion  is  fastest,  and 
energy  diffusion  appears  to  be  intermediate.  This  is  consistent  with  published  estimates  [e.g.,  Shprits 
et  a/.,  2008a, b]. 

Along  with  the  energy /spatial  patterns,  like  those  shown  in  Figures  22-24,  each  PC  has  a  time  series  of 
amplitudes.  Figures  25  and  26  show  the  time  series  of  all  eight  PC  amplitudes  along  with  observed 
conditions  for  two  of  the  so-called  GEM  storms  (named  for  NSF’s  Geospace  Environment  Modeling 
program,  which  identified  these  storms  for  coordinated  study).  In  the  October  1998  storm,  (Figure  25) 
PC-1  shows  a  gradual  increase  over  several  days  associated  with  the  main  phase  of  the  storm  (October 
~  19th),  but  preceding  the  peak  solar  wind  speed  (Vsu;).  PCs  2  and  3  make  complementary  excursions 
during  the  main  phase,  likely  associated  with  the  reconfiguration  of  the  global  magnetic  field. 

Although  the  GOES  >2  MeV  electron  flux  shows  a  fairly  straightforward  2-3  order  of  magnitude 
increase  in  (lux  superimposed  on  diurnal  variations,  the  PCs  reveal  a  far  more  complex  global  response. 
In  most  cases,  the  PCs  do  not  return  to  near  their  prestorm  values,  suggesting  global  reconfiguration  of 
the  electron  radiation  belt  topology. 

In  the  October  2000  storm  (Figure  26),  again  PC-1  shows  a  gradual  increase  over  several  days.  PCs  2 
and  3  also  show  the  same  characteristic  complementary  variations  during  the  main  phase  (October 
~  5th).  As  before,  the  PCs  reveal  a  more  complex  response  than  the  simple  increase  in  intensity  one 
would  assume  viewing  the  GOES  data. 

Statistical  analysis  of  the  PC  amplitudes  can  reveal  associations  with  other  geophysical  phenomena. 
Figure  27  shows  the  rank  order  prediction  efficiency  (the  square  of  the  rank  order  correlation 
coefficient)  between  the  posterior  PCs  and  three  other  phenomena.  The  GOES  >  2  MeV  electron  (lux  is 
correlated  with  all  of  the  first  5  PCs.  The  simplest  interpretation  of  this  result  is  that  the  environment 
at  GEO  responds  to  a  variety  of  physical  processes,  and  is  coupled  to  dynamics  at  other  locations. 
Nonetheless,  the  strong  correlation  of  GOES  with  PC  2  is  not  surprising,  given  that  PC  2  appears  to 
represent  decoupling  of  the  inner  and  outer  zones,  and  GOES  resides  far  from  the  inner  zone. 

l  he  magnetic  storm-time  disturbance  index  Dst  can  explain  up  to  about  13%  of  the  variance  of  PC-3, 
and  about  5%  of  the  variance  of  PCs  1,  2,  and  6.  1  he  simplest  interpretation  of  this  result  is  that  PCs 
1,  2,  3,  and  6  combine  to  describe  the  adiabatic  response  of  the  electron  belts  to  global  magnetic  field 
reconfigurations. 

Finally,  Figure  27  shows  that  only  PCs  1  and  3  have  much  association  with  the  solar  wind  speed. 
However,  this  analysis  accounts  only  for  the  immediate  correlation.  To  obtain  a  more  complete 
understanding  of  the  relationships  between  the  PCs  and  the  solar  wind  speed,  we  must  examine  lagged 
correlations. 
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PC-2,  a  =90°,  18.2%  of  Variance 
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Figure  23:  The  second  posterior  principal  component  of  the  TEM-2  reanalysis.  Top:  equatorially 
mirroring  electrons.  Bottom:  selected  isocontours. 
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PC-4,  aeq=90°,  10  1%  of  Variance 
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Figure  24:  The  third  and  fourth  posterior  principal  components  of  the  TEM-2  reanalysis,  for  eqnatorially 
mirroring  electrons. 
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October  1998  GEM  Storm 


Date  Starting  18-Oct-1998 

Figure  25:  Top:  TEM-2  reanalysis  posterior  principal  components  during  the  October  1998  GEM  storm. 
Lower  panels:  observed  electron  flux  at  GEO,  observed  solar  wind  speed,  and  the  magnetic  storm-time 
disturbance  index  Dst. 
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October  2000  GEM  Storm 


Figure  26:  1  EM-2  reanalysis  posterior  principal  components  and  measurements  during  the  October  2000 
GEM  storm,  in  the  format  of  Figure  25. 
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PE  of  post-PCs  with  Other  Phenomena 


Figure  27:  Prediction  efficiency  of  TEM-2  reanalysis  posterior  principal  components  with  respect  to 
GOES  >  2  MeV  electron  flux,  Dst ,  and  solar  wind  speed  (Vsw). 
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Driving  the  Posterior  PCs  with  Vsw 
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Figure  28:  Impulse  response  of  l'EM-2  reanalysis  posterior  principal  components  with  solar  wind  speed. 


Assuming  that  each  PC  has  an  independent,  possibly  time-delayed  response  to  the  solar  wind,  we  have 
computed  the  finite  impulse  response  function  for  each  PC  to  solar  wind  speed.  Figure  28  shows  how 
solar  wind  speed  drives  the  PCs.  With  the  exceptions  of  PCs  1  and  3,  all  the  PCs  have  a  strong  initial 
response  to  Of  those,  all  but  PC  4  have  a  bipolar  response  quickly  switching  sign  and 

remaining  at  the  opposite  sign  for  up  to  2  days,  ultimately  leading  to  recovery.  PCs  1  and  3  seem  to 
have  small,  but  persistently  positive  responses,  consistent  with  a  moving  average  process.  PC  4,  on  the 
other  hand,  resembles  a  classic  first-order  autoregressive  process,  with  a  quasi-exponential  impulse 
response.  Accounting  for  time  lag,  the  solar  wind  speed  describes  relatively  more  of  the  variance  of  the 
first  4  PCs  than  the  last  4,  with  PCs  1  and  4  slightly  edging  out  components  2  and  3. 

One  could  extend  the  type  of  multivariate  time  series  analysis  shown  in  Figure  28  to  develop  a  global 
radiation  belt  forecast  model,  as  opposed  to  the  GKO-specific  models  that  are  common  today  [e.g., 
Ukhorskiy  et  a/.,  2004;  Rigler  et  a/.,  2004;  Li  et  al. ,  2001;  Baker  et  a/.,  1090].  Specifically,  one  can 
attempt  to  model  or  forecast  the  8  PCs  only,  rather  than  attempting  to  model  the  thousands  of  points 
in  the  model  grid.  In  essence,  each  PC  becomes  an  index  of  radiation  belt  dynamics,  which  can  be 
modeled  in  its  own  right,  or  alongside  the  other  PCs.  A  similar  approach  is  taken  with  the  principal 
components  of  climate  variability,  such  as  the  El  Nino  Southern  Oscillation  [see,  e.g.,  Latif  et  al. ,  1998]. 

Vassiliadis  et  al.  [2002,  2003,  2005]  have  performed  a  related  analysis  on  SAMPEX/PET  data.  In 
particular,  they  have  identified  principal  regions  of  correlated  response  to  magnetospheric  arid 
interplanetary  parameters  in  one  dimension  (Lm).  I  heir  analysis  could  easily  be  applied  to  principal 
components  or  to  the  global  reanalysis  fluxes  in  3  dimensions. 


6  Summary 


We  have  developed  and  demonstrated  a  global  trapped  electron  radiation  environment  model:  TEM-2. 
This  model  comes  in  two  forms:  a  Monte  Carlo  scenario  model  and  a  Reanalysis.  These  two  elements 
prototype  the  main  components  of  future  AE9/AP9  radiation  specification  models.  Whereas  the 
quantitative  performance  of  the  models  is  often  less  than  one  would  desire  for  engineering  or  scientific 
work,  the  prototype  demonstrates  the  analysis  and  run-tirne  algorithms  that  carry  the  problem  from 
an  observed  time  series  of  unidirectional  differential  fluxes  to  an  engineering  user  application.  Further, 
the  Reanalysis,  or  “standard  solar  cycle”  provides  a  testbed  for  new  scientific  analysis  algorithms,  in 
particular  principal  component  analysis. 

We  expect  some  improvements  upon  the  algorithms  shown  here  before  the  release  of  the  AE9/AP9 
beta  and  subsequent  versions: 


•  The  AE9/AP9  beta  models  use  a  drift- invariant  A",  4>  coordinate  system  in  place  of  the  aeq ,  Lm 
system  of  TEM-2.  The  AE9/AP9  beta  will  still  use  the  OPQ  model. 

•  The  AE9/AP9  beta  will  implement  a  revision  to  equation  (20)  which  will  account  for  the 
observation  errors  explicitly  (as  opposed  to  implicitly  in  the  weighted  sum). 

•  The  AE9  beta  will  use  more  data  sets  than  were  used  in  TEM-2,  in  particular  it  will  add  LANL 
GEO  and  GPS  data. 

•  Due  to  the  finer  grid  resolution  in  AE9  and  AP9,  the  calculation  of  the  principal  components  will 
be  performed  on  a  reduced  grid  and  the  energy/spatial  patterns  interpolated  onto  a  finer  grid. 
This  extra  step  reduces  the  computational  resources  needed  to  construct  and  factor  the  spatial 
covariance  matrix  £. 

•  The  nearest  neighbors  interpolation  algorithm  will  be  modified  to  provide  an  inherently  smoother 

6. 

•  The  bootstrapping  used  to  determine  S  will  be  expanded  to  including  randomization  over 
interpolation  and  extrapolation  algorithms  and  assumptions  (e.g.,  spectral  shapes). 


The  AE9/AP9  beta,  like  TEM-2  will  still  not  have  fully- validated,  intercalibrated  data.  That,  along 
with  more  data  sets  and  the  introduction  of  a  low-altitude  grid,  will  have  to  wait  for  AE9/AP9  version 


1.0. 
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A  TEM-1 


l'EM-l  [described  briefly  in  O'Brien,  2008]  was  built  from  S3-3  MES  and  CRRES  MEA/HEEF  data 
only.  The  statistical  tables  and  scenario  generation  algorithms  were  essentially  the  same  as  those  in 
1EM-2.  However,  the  data  processing  (i.e.,  Figure  3)  was  considerably  different. 

File  nearest  neighbors  interpolation  algorithms  were  not  used.  Instead  the  flux  maps  were  built  by 
fitting  a  neural  network  to  the  observations.  Also,  0  was  defined  simply  as  (log  log  mgs),  which 
proved  to  be  a  mistake  it  did  not  explicitly  disallow  mgs  >  rriso,  especially  in  the  perturbed  maps.  In 
addition,  the  perturbations  to  0  were  problematic  with  the  neural  network  because  the  error  covariance 
matrix  for  the  network’s  free  parameters  was  often  singular,  leading  us  to  have  to  use  less  accurate 
networks  with  non-singular  error  covariance  (we  have  since  learned  how  to  circumvent  this  singularity 
by  factoring  the  error  covariance  matrix  into  its  leading  principal  components). 

The  spatial  covariance  was  obtained  by  fitting  a  simple  analytical  function,  with  an  independent 
exponential  decay  in  log  E,  aeq,  and  Lm.  We  could  not  successfully  build  neural  networks  for  the 
covariance  because  it  we  could  not  find  an  economical  way  to  constrain  either  unity  along  the  diagonal 
or  positive  definite- ness. 

In  TEM-1  we  made  some  attempts  to  intercalibrate  the  data  sets  statistically,  but  those  were 
abandoned  because  they  were  too  ad  hoc,  and  because  S3-3  MES  and  CRHES/MEA  were  already 
likely  adequately  intercalibrated  by  Vampola  on  the  ground  before  launch.  From  FEM-2  onward,  data 
intercalibration  processes  are  considered  separate  from  the  model  construction. 

A  final  lesson  that  TEM-1  taught  us  was  that  run-time  speed  depends  strongly  on  the  speed  of 
computing  the  magnetic  coordinates.  In  particular  computing  Lm  requires  tracing  a  magnetic  field 
line,  and,  looking  to  the  AE9/AP9  models,  computing  <1>  would  require  many  field  line  traces.  We 
adopted  fast  algorithms  attributed  to  K.  Pfitzer  for  computing  the  integral  invariant  /  behind  Lm 
[ Mcllwain ,  1961].  This  also  began  a  line  of  research  toward  fast  algorithms  using  analytical  functions 
and  eventually  neural  networks  [similar  to  Koller  et  al. ,  2009]. 

With  the  TEM-1  lessons  learned  in  mind,  TEM-2  was  developed  from  scratch  as  a  generic  “Next 
Generation  Radiation  Specification”,  with  algorithms  that  could  accommodate  changes  in  the  grid 
coordinates,  particle  species,  etc.  TEM-1  proved  invaluable  in  identifying  where  the  tough  problems 
were  both  from  an  algorithmic  and  a  computational  perspective. 
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B  Useful  derivatives 


In  the  numerical  optimization,  we  need  expressions  for  the  first  and  second  derivatives  of  our 
probability  equations.  We  start  with  a  restatement  of  the  negative  log  likelihood  function  from  (48) 
with  generic  marginal  distribution: 

^  -  ^2  (losx» + los/«(x«)i 

k  i 

+  ~(qJ  q  —  zr  z)  +  constants.  (92) 


The  optimization  relies  on  derivatives  of  i  with  q: 
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In  the  special  case  of  Weibull  x,  we  have: 

t  =  k)  +  \{qrq-zrz) 
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In  the  special  case  of  log-normal  x ,  we  have: 
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